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PREFACE 


Work  on  this  contract  has  examined  two  subjects,  the  analysis  of  data  using  the  NORESS  array 
in  southern  Norway  and  the  attenuation  properties  of  the  crust.  As  our  fined  report  we  present 
two  preprints  of  papers  submitted  or  in  press.  The  first,  “Phase  velocity  estimation  of  diffusely 
scattered  waves”,  by  A.  M.  Dainty  and  D.  B.  Harris,  has  been  submitted  to  the  Bulletin  of  the 
Seismological  Society  of  America.  It  deals  with  the  work  at  NORESS,  which  was  partly  carried  out 
under  a  previous  contract,  partly  in  conjunction  with  the  Lawrence  Livermore  National  Laboratory 
and  partly  under  this  contract.  The  second  paper,  “A  model  for  attenuation  and  scattering  in  the 
earth’s  crust”,  by  M.  N.  Toksoz,  A  M.  Dainty,  E.  Reiter  and  R.-S.  Wu,  is  in  press  in  PAGEOPH. 
It  describes  the  work  on  attenuation  and  was  supported  by  the  U.  S.  Geological  Survey  and  this 


contract. 
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ABSTRACT 


Two  methods  are  investigated  for  estimating  the  phase  velocity  of  diffusely  scattered  seismic 
waves  simultaneously  arriving  from  different  azimuths  and  recorded  bv  a  two-dimensional 
array  of  seismometers.  The  Hankel  spectrum  is  the  average  of  the  frequency-wavenumber 
(FK)  spectrum  over  all  azimuths,  while  the  wavenumber  spectrum  is  derived  by  integrating 
the  FK  spectrum  around  a  contour  of  constant,  phase  velocity,  i  e  .  a  circle  centered  on 
the  origin  in  the  wavenumber  plane  If  the  conventional  estimate  of  the  FK  spectrum 
using  the  covariance  matrix  of  the  seismometer  signals  is  integrated,  a  closed  form  for 
both  the  Hankel  spectrum  and  the  wavenumber  spectrum  may  be  found;  the  two  spectra 
are  very  similar,  the  wavenumber  spectrum  being  equal  to  the  Hankel  spectrum  times  the 
wavenumber  In  spite  of  this  similarity,  however,  we  find  that  the  two  formulations  have 
significantly  different  behaviour  for  small  wavenumbers,  i.e.,  high  phase  velocities  In  both 
cases  there  is  a  highest  (true)  velocity  that  can  be  estimated  from  the  spectral  maximum 
for  a  given  array  aperture  (“velocity  cutoff’’)  The  Hankel  spectrum  estimates  too  high 
a  velocity;  for  true  wavenumbers  below  a  certain  limit  infinite  velocity  is  estimated  The 
wavenumber  spectrum,  on  the  other  hand,  estimates  too  low  a  velocity,  and  there  is  an 
upper  limit  on  the  estimated  velocity.  An  example  illustrating  these  difficulties  for  the 
two  methods  is  given  for  teleseismic  P  coda  of  an  event  recorded  at  the  NORESS  array  in 
southern  Norway,  in  spite  of  the  problems,  the  analysis  is  able  to  demonstrate  that  the 


2 


coda  consists  of  two  components,  a  coherent  P  wave  component  with  a  high  phase  velocity 
and  a  diffuse  S  wave  component  of  low  phase  velocity.  The  cutoff  and  bias  problem  are 
investigated  by  numerical  simulation  for  the  NORESS  array  using  azimuthal  averaging  and 
synthetic  signals.  The  results  confirm  and  quantify  the  cutoff  problem  at  low  wavenumbers 
and  indicate  that  wavenumbers  estimated  from  the  Hankel  and  wavenumber  spectra  maxima 
bracket  the  true  wavenumber,  with  the  Hankel  spectrum  estimate  being  low  (phase  velocity 
too  high)  and  the  wavenumber  spectrum  estimate  high.  The  bias  of  both  methods  decreases 
with  increasing  wavenumber  (decreasing  phase  velocity)  and  they  are  both  asymptotically 
unbiased.  The  wavenumber  spectrum  has  a  superior  performance  at  low  wavenumbers  (high 
phase  velocity),  but  the  Hankel  spectrum  gives  superior  results  at  high  wavenumbers  (low 
phase  velocity).  The  product  of  the  linear  wavenumber  (—  ! /wavelength)  and  the  array 
aperture  define  ‘high"  and  “low"  wavenumbers;  for  low  wavenumbers  the  product  is  1 
or  less.  In  an  Appendix  we  find  absolute  lower  bounds  on  the  cutoffs  analytically.  The 
problems  could  be  mitigated  by  using  high  resolution  methods 
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INTRODUCTION 


With  a  two-dimensional  array  of  seismometers  it  is  possible  to  determine  the  phase  velocity 
(i.e.,  the  velocity  at  which  wavefronts  sweep  across  the  array)  and  azimuth  of  seismic  waves 
recorded  by  the  array  ;  Lacoss  ft  al.  (1969)  describe  some  of  the  commonly  used  techniques 
for  accomplishing  this  Generally,  most  methods  have  concentrated  on  making  as  complete 
a  determination  as  possible.  One  of  the  most  useful  concepts  has  been  estimation  of  the 
frequency-wavenumber  (FK)  spectrum,  which  is  the  power  spectrum  in  frequency  and  two 
wavenumber  components  of  the  seismic  wavefield  considered  as  a  function  of  time  t  and  two 
space  coordinates  x  and  y.  This  is  equivalent  to  decomposing  the  observed  wavefield  into 
plane  waves  with  a  range  of  azimuth  from  0  to  360"  and  a  range  of  phase  velocities  from 
0  to  oo  for  the  case  of  a  continuous  a-ray.  Specifically,  at  a  specified  angular  frequency 
ui  a  point  fcr,  ky  in  the  plane  defined  by  two  wavenumbers  corresponds  to  an  azimuth  9 
measured  from  the  y  (North)  axis  given  by 

tan  9  =■  kj  kv  ( I ) 


and  a  phase  velocity  C. 

C  —  'jj  k 

where  the  scalar  wavenumber  k  is  given  by 


(2) 


(3) 


4 


< 


Henceforth  k  will  be  referred  to  as  simply  the  wavenumber.  If  the  plane  wave  is  interpreted 
as  a  wave  of  velocity  V  incident  at  an  angle  i  to  the  vertical  (z)  axis, 

C  =  17  sin  i  (4) 

Note  that  (4)  implies 

e  >  v  (5) 

The  decomposition  into  plane  waves  represented  by  the  FK  spectrum  has  great  useful¬ 
ness  if  the  wavefield  can  be  considered  as  being  due  to  a  source  or  sources  sufficiently  far 
from  the  array  for  wavefront  curvature  to  be  ignored.  Then  significant  energy  from  azimuth 
9  can  be  interpreted  as  coming  from  a  source  situated  along  azimuth  9.  The  phase  velocity 
can  often  be  used  to  identify  the  wave  type  (from  (5))  and  calculate  its  angle  of  incidence 
for  a  body  wave  However,  some  difficulties  have  emerged  One  of  these  is  the  case  of  more 
than  one  source  or  wave  type,  especially  when  one  of  the  sources  or  wave  types  is  domi¬ 
nant  and  the  others  weak  While  some  successes  have  been  reported  (Lacoss  et  al..  1969; 
Capon  et  al.,  1969)  using  analysis  methods  that  retain  a  complete  separation  of  all  of  the 
variables  of  the  FK  spectrum,  another  possibility  is  averaging  over  some  of  the  variables  to 
increase  the  reliability  and  signal  to  noise  ratio  of  relevant  aspects  of  the  data  Nawab  et  al. 
(1985)  introduced  the  zero-delay  wavenumber  spectrum,  which  is  theoretically  equivalent 
to  summing  the  FK  spectrum  over  all  frequencies  In  addition.  Nawab  et  al.  integrated  the 
zero-delay  spectrum  along  lines  of  constant  azimutn:  by  (1)  these  are  straight  lines  passing 
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radially  outwards  from  the  wavenumber  origin.  In  this  manner  these  authors  were  able  to 
estimate  the  azimuth  of  sources  using  information  from  all  frequencies  and  phase  velocities 
present  in  the  data. 

In  this  paper  we  are  interested  in  determining  the  phase  velocity  of  waves  recorded 
at  an  array  averaged  over  all  azimuths.  This  problem  has  been  examined  by  Haubrich 
and  McCamy  (1969)  and  Iyer  and  Healy  (1972)  in  analyses  of  seismic  noise.  Haubrich 
and  McCamy  introduce  the  Hankel  spectrum,  which  is  the  average  of  the  FK  spectrum 
over  all  9  at  fixed  k  and  w.  We  have  conducted  a  theoretical  evaluation  of  the  resolution 
of  this  method  and  a  very  similar  formulation  we  have  developed  to  estimate  the  power 
in  a  window  at  time  t  as  a  function  of  wavenumber  k ;  such  an  estimate  will  be  referred 
to  as  a  wavenumber  spectrum.  The  evaluation  indicates  that  while  both  methods  have 
difficulties  at  high  phase  velocities  if  the  conventional  FK  spectrum  is  used  as  the  basis  of 
t tie  estimation,  the  wavenumber  spectrum  is  superior.  The  problem  at  high  phase  velocity 
is  due  to  the  lack  of  resolution  of  the  conventional  FK  spectrum  The  Hankel  spectrum, 
however,  gives  superior  results  for  low  phase  velocities.  We  confirm  these  results  numerically 
for  the  NORESS  array  geometry  (Bungum  et  ai.  1985;  Figure  1).  Iyer  and  Healy  (1972) 
derive  an  approximate  method  to  find  the  dominant  wavenumber  k  using  signals  from  an 
array  which  consists  of  seismometers  at  the  vertices  of  a  hexagon  surrounding  a  central 
seismometer.  We  have  not  examined  this  method  because  we  are  interested  in  a  more 
general  analysis. 
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To  illustrate  the  application  of  the  two  formalisms  we  show  an  example  using  the  seismic 


problem  of  teleseismic  P  coda.  This  example  shows  the  problems  at  high  phase  velocities 
and  is  presented  before  the  theoretical  discussion  of  these  problems  to  provide  a  practical 
demonstration  of  them.  Figure  2  shows  the  signal  received  at  the  center  seismometer  of 
the  NORESS  array  in  southern  Norway  from  a  suspected  underground  nuclear  test  in  East 
Kazahkstan  at  38°  distance.  The  large  pulse  delineated  by  the  earliest  5  s  window  is  the 
direct  P  wave  from  the  source  to  the  receiver.  A  model  for  the  energy  (“coda”)  following 
this  first  arrival  has  been  proposed  (Dainty,  1988).  The  model  hypothesizes  that  the  coda 
consists  of  two  components  due  to  scattering  in  the  crust  near  the  source  and  receiver 
respectively.  The  portion  of  the  coda  due  to  source  scattering  will  have  approximately  the 
same  phase  velocity  and  azimuth  as  the  direct  wave  and  will  be  coherent  across  the  array. 
Scattering  near  the  receiver,  on  the  other  hand,  will  produce  waves  having  the  full  range  of 
azimuths  0  to  360°  and  phase  velocities  typical  of  crustal  phases.  Since  this  component  at 
any  instant  of  time  consists  of  a  number  of  waves  of  different  azimuth  and  phase  velocity,  it 
will  not  be  coherent  across  the  array;  it  is  in  this  sense  that  it  may  be  described  as  “diffusely 
scattered”.  The  Hankel  and  wavenumber  spectra  are  applied  to  this  problem  using  the  data 
shown  in  Figure  2. 
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THEORY 


In  order  to  detect  the  presence  of  scattered  P  and  S  waves  in  the  coda,  we  propose  to 
investigate  how  coda  energy  is  distributed  in  slowness  (wavenumber)  Since  scattered  waves 
may  be  approaching  from  a  wide  range  of  azimuths,  it  is  appropriate  to  integrate  the  energy 
contributions  at  a  particular  slow'ness  over  all  azimuths.  By  scanning  over  slowness,  one 
may  attempt  to  detect  peaks  in  energy  at  the  slowness  values  which  correspond  to  P  and  S 
wave  propagation. 

Let  us  denote  the  signal  received  at  a  single  sensor  by  s(t).  Since  an  array  is  a  collection 
of  sensors,  the  received  signal  is  a  vector  waveform,  which  we  call  s(t)  We  assume  further 
that  the  data  are  available  in  digitized  form: 


sin;  =  s(tn  nAf) 


(6) 


where  we  use  the  square  brackets  to  denote  a  discrete  index,  and  At  to  denote  the  sampling 


interval. 


Following  Capon  (1969),  the  conventional  FK  spectrum  is  a  discrete  Fourier  t  ransform  in 
the  wavenumber  domain,  operating  on  the  monochromatic  spatial  rovariance  matrix  ff(w). 
with  elements  Rmn(u): 


U  M 


P{k,9) 


E 


L~.  f 


Rr, 


iki1  ( tJi  i 


"i  -  1  n  =  l 


CO 
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P(k.9)  is  defined  in  terms  of  the  azimuth  9,  the  direction  cosine  vector 

‘cost 


m  = 


sin  9 


(8) 


and  the  wavenumber  k.  The  locations  of  the  M  sensors  in  the  array  are  denoted  by  the 
vectors  i,. 

In  practice,  the  monochromatic  covariance  matrix  is  estimated  through  a  direct  Fourier 
transform  of  the  vector  signal: 


5M  =  X>, (*)£,"(*) 


(9) 


1=1 


where 


r,M  -  .4,  YL  *["  +  («  -  OA'1  e  (10) 

n=0 

It  is  common  practice  (Capon  (1969))  to  average  the  covariance  matrix  estimate  over  mul¬ 
tiple  windows  in  order  to  make  it  full-rank  and  to  reduce  the  variance  of  the  estimate.  The 
L  windows  may  be  weighted  differently  with  the  weighting  function  .4,.  Since  our  interest 
is  in  the  structure  of  the  coda,  and  since  the  amplitude  of  the  coda  decays  exponentially 
with  time  (Dainty,  1988),  we  have  made  the  weighting  inversely  proportional  to  a  decaying 
exponential  fit  to  the  envelope  of  the  coda. 


A,  -  .4  1,1  ~  DA" At 


- 1 


(11) 


We  have  chosen  o  to  fit  the  exponential  to  the  coda  envelope,  and  .4  so  that: 


12) 
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Circular  Contour  Integrals 


To  examine  the  contribution  of  diffusely  scattered  waves  to  the  coda  we  want  to  integrate 
the  FK  spectrum  (7)  around  contours  of  constant  phase  velocity.  Since  we  are  estimat¬ 
ing  a  monochromatic  FK  spectrum,  contours  of  constant  phase  velocity  are  circles  in  the 
wavenumber  plane.  Haubrich  and  McCamy  (1969)  proposed  the  Hankel  spectrum,  which 
is  the  average  of  the  FK  spectrum  over  all  9  at  a  given  k: 

H(k)=  ~  f^P(k.9)  M  (13) 

We  shall  examine  this  statistic  and  a  very  similar  one  we  shall  call  the  wavenumber  spectrum, 
which  is  proportional  to  the  energy  in  the  FK  spectrum  at  wavenumber  k  over  all  9 . 

P(fc)  =  1  P(k.  9)  kd9  (14) 

2tt  J _ w 

The  expression  kd9  in  (14)  is  the  line  element  on  the  circle  of  integration  in  the  wavenumber 
plane,  since  we  wish  to  find  the  total  energy.  Note  that,  since  k  does  not  depend  on  9 , 

P(k)  =  kH(k)  (15) 

These  integrals  can  be  evaluated  in  closed  form  for  P(k.9)  given  by  (7)  (Haubrich  and 
McCamy.  1969)  by  transforming  the  difference  in  sensor  locations  into  polar  coordinates: 

"  COs(  tl'mn  ) 

-  r„  -  Pmn  (IS) 

.Sin(Vmn)  . 
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where 

Pmn  =  \j(lm  ~  X.n)T{*m  ~  *n)  (l7) 

is  the  magnitude  of  the  difference  vector,  and  tpmn  is  the  angle  between  xm  and  xn.  Then 
using: 

dT{6)(  xm  -  xn  )  =  Pmn  cos(  0  -  )  (18) 

and  the  integral  representation  for  the  zeroth  order  Bessel  function  Jq: 

i  f'  ^1*111,-5,1  ^  =  Jo{kpmn)  (19) 

2?T  J -r 

Haubrich  and  McCamy  obtain  the  Hankel  spectrum  as. 

M  M 

H(k)  =  £  £  R  mn(u)MkP  mn  )  (20) 

m=l  n-1 

We  observe  that  (20)  can  be  computed  directly  from  the  covariance  matrix,  without  an 
intermediate  computation  of  the  FK  spectrum  itself.  Using  (15)  we  see  immediately  that 
the  wavenumber  spectrum  is  given  by. 

M  M 

P(k)  =  k  £  £  Rmn(*)Mkpmn)  (21) 

m=  1  n-  l 

Due  to  the  fact  that  the  Bessel  function  evaluations  may  be  computationally  expensive, 
one  may  wish  to  exploit  symmetries  in  the  matrix  to  reduce  the  number  of  evaluations. 
Making  use  of 

RmnM  -  K„mM  (22) 
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and 


Pnn  ~  0 


(23) 


the  number  of  evaluations  may  be  reduced  by  more  than  a  factor  of  two: 

_  M  MM 

P(k)  =  k  Y.  Rmm(v)  +  2k  £  Re{Rmn('M)}  Mkpmn)  (24) 

m=  1  m=l  n=m  1 


with  a  similar  expression  for  H(k). 
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EXAMPLE 


In  this  section  we  present  an  application  of  the  Hankel  and  wavenumber  spectra  to  the 
event  shown  in  Figure  2  and  compare  the  results  with  FK  analysis.  As  a  result  of  this 
comparison  we  note  that  the  phase  velocity  determined  by  a  finite  array  is  biased  towards 
high  values  (i.e.,  low  wavenumbers)  for  the  Hankel  spectrum  and  towards  low  values  (i.e., 
high  wavenumbers)  for  the  wavenumber  spectrum.  These  biases  effectively  constitute  a 
“velocity  estimation  cutoff”  in  the  sense  that  they  effectively  prevent  the  estimation  of  true 
velocities  higher  than  a  certain  value,  which  differs  for  the  two  spectra.  We  investigate 
the  effect  further  in  the  following  section  both  theoretically  and  using  synthetic  data  and 
numerical  techniques. 

Two  sets  of  windows,  5  s  long,  are  shown  in  Figure  2.  The  first  set  consists  of  a  single 
window  around  first  P.  The  second  set  covers  the  coda  from  20  to  100  s  after  first  P.  This 
time  interval  contains  the  arrival  PP  (surface  reflected  P),  but  the  power  in  this  phase  is  not 
large.  To  improve  the  reliability  of  the  frequency-wavenumber  and  wavenumber  spectra  for 
the  coda,  estimates  have  been  made  for  each  5  s  window  and  then  averaged  over  all  windows 
after  multiplication  of  each  estimate  by  the  weight  given  by  (II)  above.  This  procedure 
assumes  that  the  coda  is  statistically  stationary  (Dainty,  1988). 

Figure  3  shows  FK  spectra  for  first  P  and  the  coda.  In  this  Figure  and  succeeding 
Figures  and  discussion,  linear  wavenumber  K  ~  A?/(2tt)  will  often  be  used  The  spectra  are 
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taken  at  the  peak  frequency  of  the  first  P.  which  is  also  the  peak  frequency  of  the  coda.  The 
first  P  shows  a  prominent  peak  at  C  =  15.5  km/'s,  0  =  85°,  close  to  the  expected  values  of 
13.2  km/s  and  78°.  No  other  significant  energy  is  evident  in  the  first  P  window.  The  coda 
FK  spectrum  also  shows  a  peak  at  C  =  11.7  km  s  and  9  —  70°.  This  peak  indicates  there 
is  a  coherent  part  of  the  coda  due  to  scattering  near  the  source  The  difference  between 
the  phase  velocity  and  azimuth  of  the  first  P  and  the  coda  is  probably  due  to  the  combined 
effects  of  local  scattering  distorting  the  wavefrouts,  the  small  size  of  the  array,  the  use  of 
only  one  window  for  first  P.  and  contamination  of  the  coda  peak  by  energy  with  slower 
phase  velocities  and  off-source  azimuths.  There  is  some  evidence  of  a  diffuse  component  of 
the  coda  consisting  of  slower  waves  coming  from  different  azimuths,  but  it  is  not  possible  to 
determine  the  phase  velocity  or  the  relative  importance  of  this  component  from  Figure  3 
In  Figures  4  and  5  we  present  the  Hankel  and  wavenumber  spectra,  respectively,  cor¬ 
responding  to  the  FK  spectra  of  Figure  3.  Both  the  spectra  for  first  P  and  for  coda  are 
dominated  by  a  peak  at  K  =  0.  or  infinite  phase  velocity.  From  synthetic  examples,  this 
peak  is  interpreted  as  being  due  to  the  high  velocity  component  seen  in  Figure  3.  but 
the  position  of  the  peak  is  biased  towards  low  wavenumbers  (high  phase  velocities)  by  the 
weighting  effect  of  the  Bessel  function  in  (20)  towards  the  origin  of  wavenumber,  where  it 
has  a  maximum.  Careful  comparison  of  the  first  P  and  roda  Hankel  spectra  shows  that 
there  is  an  excess  of  higher  wavenumbers  in  the  coda  relative  to  first  P.  indicating  that  there 
is  some  low  phase  velocity  energy  present  in  the  coda  as  well  as  the  high  phase  velocity 
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component. 

The  wavenumber  spectra  in  Figure  5  show  a  prominent  peak  for  the  first  P  window  at 
linear  wavenumber  0.2  km'1  (C  =  9  km/s);  this  peak  is  interpreted  as  corresponding  to 
the  peak  at  C  =  15.5  km/s  seen  in  Figure  3,  but  with  a  bias  towards  lower  phase  velocities 
(higher  wavenumbers).  Note  that  due  to  the  multiplier  k  in  (15)  there  cannot  be  a  peak  at 
the  origin  of  wavenumber  for  the  wavenumber  spectrum.  There  is  a  secondary  peak  at  a 
linear  wavenumber  of  0.7  km-1  ( C  =  2.6  km/s);  this  is  apparently  a  sidelobe  of  the  main 
peak  due  to  the  finite  aperture  and  spatial  sampling  of  the  array.  This  identification  is 
confirmed  by  synthetic  examples  discussed  later.  The  coda  wavenumber  spectrum  shows  a 
peak  at  wavenumber  0.25  km-1  (C  =  7  2  km/s),  similar  to  the  values  for  the  first  P  window. 
Comparing  this  result  with  Figure  3,  the  coda  peak  at  0  25  km'1  corresponds  to  the  peak 
in  the  frequency-wavenumber  spectrum  at  C  =  11.7  km/s  and  9  =  70°,  and  represents  the 
coherent  part  of  the  coda  due  to  scattering  near  the  source.  Differences  from  the  peak  in  the 
first  P  window  are  attributed  to  the  same  factors  as  for  the  frequency-wavenumber  spectra 
of  Figure  3. 

However,  in  Figure  5  it  is  clearly  seen  that  the  coda  contains  a  component  with  higher 
wavenumbers  (lower  phase  velocities)  than  the  first  P,  with  a  broad  peak  at  a  wavenumber 
of  0.55  km-1  (C  =  3.3  km/s).  This  peak  is  not  present  in  the  first  P  window  or  for  the 
synthetic  examples  and  is  at  a  different  wavenumber  than  the  sidelobe  seen  in  the  first  P 
window  and  the  synthetic  examples.  There  is  no  single  prominent  peak  in  Figure  3  with  this 
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phase  velocity.  A  phase  velocity  of  3.3  km/s  is  characteristic  of  horizontally  propagating  S 
waves  or  surface  waves.  To  arrive  so  soon  after  first  P,  such  waves  must  have  been  generated 
near  the  array.  This  peak  is  thus  identified  with  a  diffusely  scattered  component  of  the  coda 
consisting  of  shear  or  surface  waves  arriving  from  many  different  azimuths  due  to  scattering 
with  phase  velocities  between  3  and  4  km/s.  Since  the  wavenumber  spectrum  gives  an 
estimate  of  the  energy  present  at  a  given  wavenumber,  comparing  first  P  and  coda  spectra 
in  Figure  5  indicates  that  the  diffusely  scattered  component  appears  to  be  about  as  strong 
as  the  coherent  component  for  this  event. 

In  summary,  the  analysis  using  the  wavenumber  spectrum  has  successfully  separated  the 
coda  into  near  source  and  near  receiver  scattered  components  However,  a  bias  towards  high 
wavenumbers  (low  phase  velocity)  was  noted.  The  FK  spectrum  analysis  did  not  separate 
the  components  clearly  because  of  the  diffuse  nature  of  the  near  receiver  component  The 
Hankel  spectrum  analysis  encountered  difficulty  due  to  a  bias  towards  low  wavenumbers 
(high  phase  velocity)  which  led  to  domination  of  the  spectra  by  the  high  phase  velocity  near 
source  component.  In  the  next  section  we  examine  the  issue  of  bias,  which  first  became 
evident  in  analyses  such  as  those  presented  above 
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VELOCITY  ESTIMATOR  BIAS 


Both  the  Hankei  and  wavenumber  spectra  can  be  used  as  estimators  of  phase  velocity.  Since 
they  are  estimators  of  average  or  total  energy  respectively  as  a  function  of  wavenumber,  by 
finding  a  maximum  of  energy,  we  obtain  an  estimate  of  wavenumber  (and  for  monochromatic 
signals,  phase  velocity).  The  performance  of  estimators  is  usually  characterized  by  two 
quantities:  bias  and  variance.  In  this  section  we  address  the  question  of  bias,  leaving  the 
variance  to  a  future  publication 

The  Hankei  and  wavenumber  integrals  have  nearly  identical  form,  but  exhibit  different 
bias.  We  use  both  because  the  estimates  they  provide  bracket  the  true  wavenumbers  (and 
phase  velocities)  of  incident  waves.  In  addition,  we  demonstrate  that  bias  is  controlled 
by  array  aperture.  For  arrays  of  limited  aperture,  the  bias  problem  becomes  acute  at 
small  wavenumbers,  culminating  in  an  estimation  cutoff  The  existence  of  a  cutoff  means 
that  wavenumbers  below  a  calculable  threshold  cannot  be  estimated.  In  an  Appendix, 
an  absolute  lower  bound  on  this  threshold  which  depends  only  on  aperture  is  obtained 
analytically  for  arrays  exhibiting  near-circular  symmetry.  In  this  section  we  numerically 
evaluate  the  expected  bias  and  the  exact  cutoff  for  the  NORESS  array  geometry,  and 
confirm  the  calculations  using  synthetic  signals. 

The  cutoff  phenomenon  is  examined  by  considering  the  spectra  generated  by  plane  waves 
incident  on  the  array.  Suppose  a  plane  wave  is  incident  from  direction  0.,  with  angular 
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frequency  w  and  phase  velocity  C .  The  theoretical  covariance  matrix  for  the  wavefield  is: 

«mn(w)  =  E  e  ^  (25) 

where  E  is  the  energy  density  of  the  wave.  The  Hankel  and  wavenumber  spectra  are: 

P(k)  =  kE  ]P  e  C-  o)  -m  -n)  Jn{kpmn)  kH(k)  (26) 

m=  1  r»  =  l 

These  spectra  should  be  approximately  independent  of  the  direction  of  approach.  90. 
for  arrays  with  nearly  circularly  symmetric  geometry  such  as  NORESS  (Figure  1).  Indeed. 
NORESS  was  designed  to  have  uniform  performance  over  all  azimuths  Instead  of  picking 
a  particular  direction  for  our  analysis,  we  examine  the  average  spectra: 


1  r  - 

2  ,/_/<*> 


k  H{k) 


Because  of  the  rotational  symmetry  of  the  problem,  the  average  spectra  should  be  approx¬ 
imately  equal  to  the  spectra  obtained  for  a  particular  direction  of  approach:  consequently, 
the  estimated  wavenumbers  and  biases  should  be  approximately  the  same  The  average 


spectra  are: 


st  st 


**  £  E  Mb  Pmrx)  rn  n  ) 


k  H(k)' 


where 
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From  these  analytic  expressions  for  the  spectra,  the  wavenumbers  that  would  be  esti¬ 
mated  from  the  spectral  maxima  for  incident  plane  waves  having  true  wavenumber  k\  or 
true  linear  wavenumber  K' ,  can  be  calculated  numerically  for  the  NORESS  geometry.  In 
Figure  6,  these  maxima  wavenumbers  for  the  azimuthally  averaged  Hankel  and  wavenum¬ 
ber  spectra  are  plotted  as  a  function  of  K *  over  the  range  0  to  1  cycle/km.  The  cutoff 
effect  discussed  earlier  is  clearly  seen  near  K*  —  0:  for  the  Hankel  spectrum  a  value  of  0 
is  estimated  for  any  value  of  K *  less  than  ~  0.22,  while  for  the  wavenumber  spectrum  the 
estimated  value  at  K'  —  0  is  ~  0.17.  Indeed,  generally  the  wavenumber  spectrum  estimate 
is  biased  high  and  the  Hankel  spectrum  estimate  is  biased  low.  As  previously  noted,  they 
jointly  bracket  the  true  wavenumber. 

To  assure  that  the  azimuthally  averaged  spectra  provide  wavenumber  estimates  close 
to  those  obtained  for  individual  plane  waves,  we  computed  estimates  for  synthetic  signals 
representing  single  plane  waves  with  a  dominant  frequency  of  1  or  2  Hz  incident  from 
due  North  with  a  specified  phase  velocity.  Figure  7  shows  the  Hankel  and  wavenumber 
spectra  for  a  2  Hz  signal  with  a  phase  velocity  of  23  km/s.  Note  the  similarity  of  the 
results  for  this  high  phase  velocity  signal  to  the  first  P  spectra  in  Figures  4  and  5.  Maxima 
estimates  for  a  suite  of  phase  velocities  from  such  synthetic  cases  are  plotted  on  Figure 
6  and  agree  well  with  the  averaged  spectra  results,  certifying  the  averaging  approach.  In 
fact,  the  approach  of  using  the  azimuthally  averaged  spectra  may  be  more  correct  for  coda 
scattering  investigations,  since  waves  may  be  expected  to  arrive  over  a  band  of  directions. 
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The  implications  of  this  analysis  for  velocity  estimation  are  more  readily  appreciated  in 
Figure  8.  In  this  Figure,  the  wavenumber  estimates  are  reinterpreted  as  velocity  estimates 
for  the  two  cases  /  =  1.0  Hz  and  /  -  2.0  Hz.  The  true  velocity  is  bounded  from  above 
by  the  Hankel  spectrum  estimate  and  from  below  by  the  wavenumber  spectrum  estimate. 
Low  velocities  are  estimated  well,  but  the  spectra  provide  no  constraints  in  the  high  ve¬ 
locity  range.  With  increasing  frequency  the  constraints  improve,  i  e..  velocity  resolution 
increases  with  frequency,  because  at  a  given  phase  velocity  the  wavenumber  is  higher  at 
higher  frequencies  by  (2) 

As  mentioned  previously,  the  cutoff  problem  is  examined  analytically  in  an  Appendix. 
By  setting  pmn  =  A.  where  A  is  the  diameter  (aperture)  of  the  array,  it  may  be  shown  that 
the  Hankel  spectrum  will  have  a  peak  at  k  -  0  for  all 

•-  o  A  (32) 

where  a  =s  2.4.  Since  actually  pmn  <  A.  this  is  a  lower  bound  on  the  cutoff  for  the  Hankel 
spectrum.  Similarly,  the  wavenumber  spectrum  cannot  have  a  peak  at  wavenumbers  below 

*mm  d  A  (33) 

where  d  1.2.  Again,  this  is  a  lower  bound.  Remembering  that  k  jj  C .  the  maximum 
phase  velocity  that  can  be  estimated  by  the  wavenumber  spectrum  is  determined  from: 

^A  •  ^  (34) 
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For  the  Hankel  spectrum,  the  cutoff  is  defined  by: 
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max 


u>  A 
2~4 


(35) 


(36) 


The  two  theoretical  wavenumber  cutoffs  implied  by  (32)  and  (33)  are  indicated  on  Figure 
6.  It  is  clear  that  the  actual  NORESS  cutoffs  occur  at  higher  wavenumbers.  This  confirms 
that  the  minimum  theoretical  wavenumbers  are  lower  bounds  which  will  not  be  attained 
by  any  realizable  array.  The  disparity  occurs  because  in  a  realizable  array  all  of  the  pmn 
cannot  be  equal  to  the  array  aperture  as  was  assumed  in  the  analysis.  However,  although 
the  bounds  are  optimistic,  they  do  establish  the  existence  of  the  cutoff  phenomenon. 

The  cutoffs  can  be  understood  intuitively  with  the  help  of  Figure  9.  This  Figure  shows 
the  FK  spectrum  for  a  pure  plane  wave  coming  from  the  North  with  a  true  linear  wavenum¬ 
ber  K*  =  0.1  cycles/km,  within  both  cutoffs.  The  spec  ,  am  is  estimated  by  computing 
the  beampattern  and  shifting  it  by  0.1  cycles/km  Northward  to  emphasise  the  influence  of 
the  width  of  the  main  lobe  of  the  beampattern  on  the  cutoff.  The  dashed  circles  represent 
contours  of  integration  in  (13)  and  (14).  The  inner  circle  passes  through  the  peak  of  the 
FK  spectrum;  if  there  were  no  cutoffs  and  no  bias,  both  the  Hankel  spectrum  (13)  and  the 
wavenumber  spectrum  (14)  would  have  their  maxima  for  this  circle  of  integration.  The  outer 
dashed  circle  is  at  the  cutoff  radius  for  the  wavenumber  spectrum,  K  =  0.17.  Considering 


2 


first  the  Hankel  spectrum,  the  value  obtained  for  the  inner,  true,  circle  will  be  the  average 
around  the  circle  by  (13).  But  we  see  that  this  is  probably  less  than  the  value  at  the  origin, 
K  —  0,  because  the  value  at  the  origin  is  ~  500,  and  most  of  the  circle  passes  through 
regions  of  the  FK  spectrum  that  have  lower  values  (the  peak  of  the  spectrum  is  normalized 
to  the  square  of  the  number  of  array  elements,  252  =  625).  Our  analysis  using  (30)  and 
synthetic  signals  presented  in  Figure  6  shows  that  in  fact  the  peak  of  the  Hankel  spectrum 
lies  at  the  origin.  For  the  wavenumber  spectrum,  the  integration  around  the  dashed  circular 
contours  is  weighted  by  the  factor  k,  the  radius  of  the  circle  of  integration.  Because  of  this, 
it  is  quite  possible  that  the  largest  value  of  the  wavenumber  spectrum  will  lie  at  a  greater 
value  of  K  than  the  true  value  shown  by  the  inner  dashed  circle  in  Figure  9.  The  results 
summarized  in  Figure  6  show  that  the  largest  value  actually  occurs  for  the  outer  dashed 
circle. 

The  crucial  factor  in  the  argument  presented  in  the  preceding  paragraph  is  the  role 
played  by  the  finite  width  of  the  main  lobe  of  the  beampattern  caused  by  the  finite  aperture 
of  the  array.  If  the  array  had  infinite  aperture,  the  broad  peak  seen  in  Figure  9  would  become 
a  very  narrow  peak  concentrated  near  the  true  value  indicated  by  the  cross.  Then  both 
the  Hankel  and  the  wavenumber  spectra  would  have  their  peaks  at  the  true  wavenumber 
because  there  would  only  be  significant  power  near  this  wavenumber.  Instead,  the  lack  of 
resolution  of  the  actual  array  causes  energy  in  the  FK  spectrum  to  be  spread  to  lower  and 
higher  wavenumbers,  with  the  results  noted  above  for  sufficiently  small  wavenumbers.  We 
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shall  return  to  this  question  of  resolution  later. 

One  final  method  of  viewing  the  bias  problem  is  instructive.  In  Figure  10  we  have 
plotted  the  normalized  velocity  bias,  defined  as: 

Ceft  -  Ciruc  _  1  Ktrt  ~  I  K' 

Ctru'  "  I  K-  (  ’ 

as  a  function  of  the  aperture  expressed  as  a  fraction  of  horizontal  wavelength  A,  A/ A  =  AK' . 
The  Hankel  spectrum  bias  exhibits  the  cutoff  very  clearly  with  a  vertical  asymptote  at  small 
aperture.  The  biases  of  both  methods  become  acceptably  small  as  the  aperture  exceeds  one 
wavelength,  which  is  a  restatement  of  the  classical  Rayleigh  resolution  limit.  Below  one 
wavelength,  the  wavenumber  spectrum  has  smaller  bias,  and  may  be  usable  down  to  about 
half  a  wavelength.  The  relative  performance  of  the  Hankel  and  wavenumber  spectra  at  low 
wavenumbers  agrees  with  the  results  (35)  and  (36)  as  well  as  our  practical  experience  as 
shown  in  the  Example  section.  Above  one  wavelength,  the  Hankel  spectrum  generally  has 
less  bias,  though  both  methods  are  asymptotically  unbiased. 

A  final  important  point  that  must  be  made  is  that  the  analysis  of  this  velocity  cutoff 
has  been  explicitly  carried  out  for  the  formulations  (24)  and  (20)  which  were  based  on  the 
conventional  FK  estimate  (7),  and  would  not  apply  to  the  contour  integrals  (13)  and  (14) 
based  on  some  other  FK  estimate.  Indeed,  both  methods  are  constrained  by  the  classical 
Rayleigh  resolution  limit  The  existence  of  higher  resolution  methods,  such  as  the  Maximum 
Likelihood  Method  (Capon,  1969),  suggests  that  the  constraint  may  be  eased,  reducing  the 
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bias,  and  extending  the  analysis  to  higher  velociti 
future  publication. 
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DISCUSSION  AND  CONCLUSIONS 


Two  methods  have  been  evaluated  to  find  the  phase  velocity  of  waves  from  two-dimensional 
array  data  when  the  wavefield  consists  of  waves  arriving  simultaneously  from  different  az¬ 
imuths.  The  first  method,  the  Hankel  spectrum  of  Haubrich  and  McCamy  (1969),  is  the 
average  of  the  FK  spectrum  over  all  azimuths  at  a  constant  wavenumber.  The  second 
formulation,  called  the  wavenumber  spectrum,  is  derived  by  integrating  the  FK  spectrum 
around  circles  centered  on  the  origin  in  wavenumber  space,  i.e.,  contours  of  constant  phase 
velocity.  If  the  conventional  estimate  of  the  FK  spectrum  using  the  covariance  matrix  of 
the  seismometer  signals  is  used,  the  resulting  integrals  may  be  evaluated  in  closed  form. 
The  two  expressions  that  result  are  very  similar:  if  k  is  the  wavenumber,  the  wavenumber 
spectrum  is  k  times  the  Hankel  spectrum 

The  spectra  have  been  used  to  examine  teleseismic  P  coda  using  the  NORESS  array. 
The  wavenumber  spectrum  demonstrates  that  the  coda  consists  of  a  high  phase  velocity 
P  wave  component  scattered  near  the  source  and  a  low  phase  velocity  S  wave  component 
scattered  near  the  receiver.  However,  a  bias  effect  at  high  phase  velocities  was  noted;  for 
the  Hankel  spectrum  there  was  also  a  bias  of  such  severity  that  useful  results  could  not  be 
obtained.  The  sense  of  the  bias  was  different  for  the  two  formulations,  leading  to  high  phase 
velocities  for  the  Hankel  spectrum  and  low  phase  velocities  for  the  wavenumber  spectrum 
For  both  formulations  the  problem  appeared  to  culminate  in  a  velocity  cutoff,  i.e..  there 
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was  a  highest  (true)  phase  velocity  that  could  be  estimated 

To  investigate  the  problem  further  and  explain  the  different  results  obtained  by  the  two 
methods  we  undertook  theoretical  and  numerical  studies.  The  apparently  slight  difference 
in  the  formulations  has  a  pronounced  effect  near  k  =  0.  i.e.,  high  phase  velocities  The 
wavenumber  spectrum  is  zero  at  k  -  0  due  to  the  factor  of  k  noted  above,  but  the  Hankel 
spectrum  has  a  tendency  to  have  a  spurious  peak  at  k  —  0  due  to  the  weighting  of  the 
zeroth  order  Bessel  function  that  appears  in  the  formulation  This  effect  was  demonstrated 
numerically  by  considering  the  azimuth  averaged  spectra  and  simulations  using  synthetic 
signals,  and  investigated  theoretically  in  terms  of  the  highest  phase  velocity  that  can  be 
determined,  i.e.,  the  velocity  cutoff.  For  the  Hankel  spectrum,  the  problem  was  to  find 
the  minimum  k  at  which  a  peak  can  be  distinguished  from  a  peak  at  k  -  0.  while  for 
the  wavenumber  spectrum  the  minimum  value  of  k  at  which  a  peak  can  occur  is  required. 
The  numerical  results  confirmed  and  quantified  the  existence  of  the  bias  and  the  velocity 
cutoff  for  the  NORESS  geometry.  The  theoretical  investigation  demonstrated  that  for 
both  formulations  there  is  a  limit  that  depends  only  on  the  array  aperture  for  circular 
geometries  The  numerical  and  theoretical  investigations  showed  that  the  wavenumber 
spectrum  is  superior  to  the  Hankel  spectrum  for  distinguishing  high  phase  velocities  but 
the  Hankel  spectrum  is  more  accurate  at  low  phase  velocities.  Both  formulations,  however, 
are  asymptotically  unbiased  for  the  limit  of  low  phase  velocities. 

It  should  be  noted  that  the  cutoff  is  due  to  the  low  resolution  of  the  conventional 
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FK  formulation  for  arrays  of  limited  aperture.  To  ameliorate  or  eliminate  this  problem 
a  number  of  solutions  could  be  attempted.  For  waves  coming  from  a  single  azimuth,  the 
simplest  solution  is  to  use  the  original  FK  formulation,  which  exhibits  no  biases,  instead  of 
the  wavenumber  spectrum.  However,  the  same  lack  of  resolution  has  limited  the  usefulness 
of  FK  when  more  than  one  wavetype  or  azimuth  is  present  Another  possibility  would 
seem  to  be  an  increase  in  array  aperture,  but  in  practice  the  prospects  for  seismic  waves 
are  considerably  limited  by  loss  of  coherence  of  the  signal  over  large  intersensor  spacings 
(Ingate  et  al.,  1985).  A  much  more  promising  solution  would  seem  to  be  to  start  from  a 
high-resolution  estimate  of  the  FK  spectrum,  such  as  the  Maximum  Likelihood  Method 
(MLM;  Capon,  1969)  MLM  has  traditionally  been  applied  to  seismic  data  with  great 
success  In  addition  to  ameliorating  the  cutoff  problem  by  its  high  resolution,  MLM  will 
also  help  the  sidelobe  problem  briefly  alluded  to  in  the  Example  but  not  further  discussed. 
One  disadvantage  of  starting  from  a  MLM  formulation  of  the  FK  spectrum  is  that  the 
contour  integrals  (13)  and  (14)  must  be  carried  out  numerically.  We  will  report  on  the  use 
of  MLM  in  a  later  publication. 
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APPENDIX 


In  this  Appendix  we  examine  the  question  of  estimation  cutoffs  analytically  starting  from 
(30).  For  the  Hankel  spectrum  the  numerical  calculations  presented  in  the  text  seemed  to 
show  that  below  a  certain  limit  (cutoff)  in  the  true  wavenumber  K'  (or  k')  the  Hankel 
spectrum  always  seemed  to  indicate  a  wavenumber  of  0.  We  argue  the  existence  of  the  limit 
by  showing  analytically  that  for  sufficiently  small  k' .  \H(k ))  is  always  maximum  at  the 
origin.  Now: 

MM 

( H(k ))  =  E  £  Y,  Mk'Pmn)  MkPmn)  (A-l) 

m=  I  n= 1 
M  M 

<  e  y.  E  Pmn  )l  \Jn{kp  mn  )  I  (A-2) 

m=l  n=l 

Jo{x)  is  maximum  at  the  origin  and  is  a  monotonically  decreasing  function  in  the  interval 
[0.  o|,  where  a  =;  2.40  is  the  first  zero  of  Jn(x).  The  intersensor  distances  pmn  are  bounded 
by  the  array  aperture  A: 

Pmn  1  A  (A  -  3) 

If  we  choose 

k'  A  <  a  (A  -  4) 

then 

■J'A't'Pmn)  ./n(fr'A)  ->  0.  (A -5) 
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Thus,  we  have: 


M  M 

(Tim  <  SEE  ^o(jb*Pmn)  |Jo(A:Pmn)|  (A-6) 

m=l  n=l 

<  £EE  M^pmn)  =  <F(0))  (A-7) 

m=l  rv=  1 

which  is  the  desired  result.  This  result  demonstrates  that  the  cutoff  phenomenon  is  depen¬ 
dent  on  aperture:  for  any  it*  <  a/ A,  the  estimated  wavenumber  will  be  0.  The  Hankel 
spectrum  cannot  distinguish  small  non-zero  wavenumbers  from  zero  wavenumber. 

The  cutoff  phenomenon  is  slightly  more  complicated  to  establish  for  the  wavenumber 
spectrum.  In  this  case  we  show  that  the  spectrum  cannot  have  a  peak  below  a  certain  cutoff 
wavenumber,  kmin.  We  establish  that: 

-j-  (P{k))  >0  for  k  <  kmin  (A  -  8) 

Now: 

,  MM 

—  (P(k))  =  E  Jo(k*Pmn)  ( MkPmn )  "  kPmn  Jl  (kPmn))  (A  -  9) 

aK  m=l  n=l 

Again,  the  first  zero  of  Jq(x )  occurs  at  a,  and  the  first  zero  of  ~  -r-M1)  occurs  at 
0  ss  1.2.  For  i  smaller  than  the  appropriate  value,  both  terms  are  positive.  If  we  choose 

it*  A  <  a  (A -10) 

then  the  derivative  is  positive  for  all  k  satisfying 

kpmn  <  kA  <  0  (A -11) 
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which  indicates  no  peaks  for  k  in  the  range  jO./?/A)- 
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FIGURE  CAPTIONS 


Figure  1.  Map  of  NORESS  array  seismometer  sites.  Array  center  is  at  60.7353°N,  11.5414°E. 

Figure  2.  Teleseismic  event  recorded  at  the  center  seismometer,  NORESS.  Solid  vertical  lines 
delineate  analysis  window  around  first  P,  dashed  lines  indicate  coda  windows.  Event  is  a 
suspected  nuclear  test  at  the  East  Kazakh  test  site  on  February  10  1985. 

Figure  3.  FK  spectra  of  the  event  shown  in  Figure  2  at  a  frequency  of  1.6  Hz.  (Top)  Five 
second  window  around  first  P.  (Bottom)  Window  between  20  and  100  seconds  after  first  P 
(coda),  calculated  as  average  of  non-overlapping  5  second  windows.  From  Dainty  (1988). 

Figure  4.  Hankel  spectra  of  the  event  shown  in  Figure  2  at  a  frequency  of  1.8  Hz.  Windows 
as  for  Figure  3. 

Figure  5.  Wavenumber  spectra  of  the  event  shown  in  Figure  2  at  a  frequency  of  1.8  Hz. 
Windows  as  for  Figure  3 

Figure  6.  Linear  wavenumber  estimated  from  average  Hankel  and  wavenumber  spectra  maxima 
by  (30)  plotted  against  true  linear  wavenumber  (solid  lines).  Long  dashed  line  is  Estimated 
Wavenumber  =•  True  Wavenumber.  Symbols  are  results  from  simulations  such  as  those 
shown  in  Figure  7  inverted  triangles  are  for  a  dominant  frequency  of  1  Hz  and  crosses  are 
for  a  dominant,  frequency  of  2  Hz.  Arrows  show  positions  of  theoretical  wavenumber  limits 
implied  by  (32)  for  the  Hankel  spectrum  on  the  true  wavenumber  axis,  and  by  (33)  for  the 
wavenumber  spectrum  on  the  estimated  wavenumber  axis 
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Figure  7.  Ilankel  and  wavenumber  spectra  for  a  synthetic  plane  wave  signal  of  phase  velocity 
23  km/s,  azimuth  due  North  and  dominant  frequency  2  Hz  placed  on  the  sensors  of  the 
NORESS  array.  Ordinate  is  linear  wavenumber. 

Figure  8.  Results  shown  in  Figure  6  replotted  as  estimated  velocity  against  true  velocity  for 
frequencies  of  1  Hz  (solid  lines)  and  2  Hz  (short  dash  lines).  Long  dash  line  is  Estimated 
Velocity  =  True  Velocity.  Symbols  as  for  Figure  6. 

Figure  9.  FK  spectrum  for  a  pure  plane  wave  incident  on  the  NORESS  array  from  the  North 
with  a  linear  wavenumber  of  0.1  cycles 'km.  Dashed  circles  represent  contours  of  integration 
in  (13)  and  (14).  See  text  for  discussion. 

Figure  10.  Results  of  Figure  6  replotted  as  fractional  velocity  bias  (37)  against  K'  A  =  A/A. 


Symbols  as  for  Figure  6. 
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Abstract 


The  mechanisms  contributing  to  the  attenuation  of  earthquake  ground  motion  in  the  distance  range 
of  10  to  200  km  are  studied  with  the  aid  of  laboratory  data,  coda  waves,  Rg  attenuation,  strong 
motion  attenuation  measurements  in  the  northeast  United  States  and  Canada,  and  theoretical 
models.  The  frequency  range  1-10  Hz  has  been  studied.  The  relative  contributions  to  attenuation 
of  anelasticity  of  crustal  rocks  (constant  Q),  fluid  flow  and  scattering  are  evaluated.  Scattering  is 
found  to  be  strong  with  an  albedo  Bo  =  0.8-0. 9  and  a  scattering  extinction  length  of  17-32  km. 
The  albedo  is  defined  as  the  ratio  of  the  total  extinction  length  to  the  scattering  extinction  length. 
The  Rg  results  indicate  that  Q  increases  with  depth  in  the  upper  kilometer  or  two  of  the  crust,  at 
least  in  New  England.  Coda  Q  appears  to  be  equivalent  to  intrinsic  (anelastic)  Q  and  indicates 
that  this  Q  increases  with  frequency  as  Q  =  Qofn,  where  n  is  in  the  range  0.2-0.9.  The  intrinsic 
attenuation  in  the  crust  can  be  explained  by  a  high  constant  Q  (500  <  Qo  <  2000)  and  a  frequency 
dependent  mechanism  most  likely  due  to  fluid  effects  in  rocks  and  cracks.  A  fluid-flow  attenuation 
model  gives  a  frequency  dependence  (Q  ~  Qof° 5)  similar  to  those  determined  from  the  analysis  of 
coda  waves  of  regional  seismograms.  Q  is  low  near  the  surface  and  high  in  the  body  of  the  crust. 
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Introduction 


The  dependence  of  ground  motion  U  of  frequency  w  with  distance  r  from  the  source  has  been 
described  by  the  equation 

U(r,  w)  =  S(w)  ■  F(r)  •  exp[-ro,/(2 Qv)}  (1) 

In  (1)  S  represents  the  source  and  u  is  the  seismic  wave  velocity.  The  factor  F(r)  is  the  geometrical 
spreading  factor  due  to  the  propagation  of  motion  in  a  perfectly  elastic  medium  without  random 
scattering.  In  a  homogeneous  medium  F[r)  =  1/r  for  body  waves  and  =  1/%/r  for  surface  waves 
(not  including  the  effect  of  dispersion),  and  this  form  is  often  used  at  ranges  of  less  than  100 
km.  The  exponential  term  parameterised  by  Q  describes  the  effect  of  attenuation  due  to  imperfect 
elasticity  (anelastic  attenuation)  and  random  scattering.  However,  a  recently  developed  theory  of 
wave  scattering  in  absorbing  and  heterogeneous  media  (Wu  1984,  1985)  indicates  that  the  spectral 
intensity  of  ground  motion  varies  with  distance  in  a  much  more  complicated  manner  when  multiple 
scattering  is  considered.  In  this  theory,  the  shape  of  the  intensity-distance  curve  is  controlled 
by  the  seismic  albedo  Bo  of  the  medium  discussed  below,  as  shown  in  Figure  1.  Only  for  the 
case  of  Bq  <?C  0.5  can  the  curve  be  approximately  expressed  in  a  simple  exponential  form  after 
geometric  spreading  corrections.  Our  objective  is  to  examine  the  relative  contribution  of  scattering 
and  anelastic  attenuation  in  the  crust  and  to  propose  a  general  crustal  model  of  attenuation  and 
scattering. 

A  hile  physical  definitions  of  Q  are  given  in  the  literature  in  practice  Q  is  usually  measured 
by  methods  based  on  equation  (1)  or  the  corresponding  time  equation.  A  difficulty  that  arises 
with  this  ‘operational’  definition  of  Q  is  that  different  results  may  be  obtained  based  on  different 
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methods  of  measurement  even  in  simple  cases  (Aki  and  Richards,  1980,  p.  168).  A  particularly 
important  case  is  the  influence  of  scattering.  Without  scattering,  the  Q  measured  using  (1)  will 
represent  solely  anelastic  attenuation.  However,  if  a  short  pulse  travels  through  a  medium  in  which 
both  scattering  and  anelasticity  occur,  then  the  total  Q,  Qt,  is  (Dainty,  1981) 

i/Qt  =  i/Q*  +  i/Q,  (2) 

In  (2)  Qa  is  the  anelastic  Q.  The  loss  of  energy  from  the  pulse  due  to  scattering  is  described  by 
Qs  —  oj/(2rjsv),  where  rj,  is  the  scattering  coefficient  discussed  below.  Equation  (2)  indicates  that 
for  the  case  of  a  short  pulse  the  influence  of  anelastic  attenuation  and  scattering  cannot  be  separated 
unless  some  model  for  the  Q  components  is  used.  However,  the  scattered  energy  represented  by  Q, 
is  merely  deflected,  not  lost  from  the  wavefield.  Measurements  of  Q  which  examine  a  substantial 
window  in  time  of  the  seismogram,  such  as  coda  Q,  may  include  some  of  this  energy. 

Below  we  present  summaries  of  three  investigations  that  measure  Q  in  different  ways.  The 
first  interprets  the  attenuation  of  strong  motion  measurements  using  pseudo  velocity.  Since  an 
integration  over  the  duration  of  strong  motion  is  carried  out  in  this  measure,  the  total  energy  in 
the  seismogram  including  scattered  energy  is  measured.  Our  investigation  explicitly  examines  the 
question  of  the  relative  attenuation  due  to  anelasticity  and  scattering  and  separates  the  contribution 
of  each  using  transport  theory.  The  second  investigation  examines  the  Rg  phase  (fundamental 
mode  Rayleigh  wave)  in  the  frequency  range  1-10  Hz  at  distances  of  10-30  km.  This  is  a  pulse 
measurement  and  equation  (2)  should  apply,  i.  e.,  the  total  Q  represents  the  combined  effect  of 
anelasticity  and  scattering  without  separating  them.  Finally  we  discuss  the  implications  of  some 
previous  measurements  of  coda  Q.  While  this  method  has  become  a  standard  method  of  measuring 
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Q,  the  interpretation  of  the  results  has  not  been  clear. 


We  consider  primarily  attenuation  in  the  distance  range  of  10  to  100  km.  This  interval  is  ideal 
for  several  important  reasons.  At  distances  shorter  than  about  10  km  from  the  source,  non-hnear 
behavior  of  materials  due  to  high  strains  ( e  >  10-5)  can  dominate.  At  distances  greater  than 
100  km,  the  geometric  spreading  effects,  due  to  velocity-depth  functions  and  multiple  branches 
of  travel-time  curves,  become  site-specific  and  uncertain.  Another  important  factor  for  favoring 
this  distance  range  is  that  a  considerable  amount  of  new  attenuation  data  has  been  obtained  both 
from  strong  motion  records  and  the  analysis  of  seismic  coda  waves.  Also,  because  of  the  rapid 
attenuation  of  high  frequency  rig  it  is  only  easily  observed  in  this  distance  range. 

Before  reviewing  the  attenuation  data,  it  is  important  to  define  the  terminology.  The  attenua¬ 
tion  for  a  given  wave  type  (P  or  S)  is  defined  as  the  inverse  of  the  quality  factor  Q,  and  related  to 
other  measures  by: 

1  _  at>  _  8 

Q  * 

where  a  is  the  anelastic  attenuation  coefficient,  v  the  wave  velocity,  /  the  frequency,  and  6  the 
logarithmic  decrement. 

Attenuation  Q~l  or  the  quality  factor  Q  are  dimensionless  quantities.  Physically, Q-1  is  equal  to 
the  ratio  of  energy  dissipated  per  cycle  to  the  total  energy.  For  small  attenuation,  (i.e.  Q~l  <  0.1), 
additional  relationships  can  be  established  in  terms  of  stress-strain  relationships: 


1 

Q 


——  =  Un<p  ~  4> 
Mr 


(4) 


where  M[  and  Mr  are  the  imaginary  and  real  parts  of  the  appropriate  elastic  modulus  (M  = 
Mr  +  iMi)  and  is  the  phase  lag  of  the  strain  behind  the  stress  (i.e.,  loss  tangent).  The  dimension 
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of  the  anelastic  attenuation  coefficient  a  is  generally  given  as  dB/unit  length  or  nepers  per  unit 
length.  The  relationship  between  the  two  is  a(dB/unit  length)  =  8.686  a(nepers/unit  length). 

Most  of  the  data  for  crustal  attenuation  comes  from  coda  waves  (Aki  and  Chouet,  1975;  Aki, 
1980;  Pulli,  1984;  Singh  and  Herrmann,  1983;  Singh,  1985;  Gupta  et  al.,  1983;  Rautian  and  Khal- 
turin,  1978;  Roecker  et  al.,  1982;  Herrmann,  1980).  These  measurements  generally  give  attenuation 
that  decreases  with  frequency  in  the  frequency  range  of  /  =0.5  to  25  Hz.  A  typical  result  for  New 
England  is  Qc{f)  —  460 f0A  (Pulli,  1984).  The  increase  of  Q  with  frequency  and  the  high  values 
( Q  >  1000)  at  frequencies  above  10  Hz  in  the  Eastern  United  States  cannot  be  reconciled  with  the 
laboratory  measurements  of  Q  in  crustal  rocks  (see  Toksoz  and  Johnston,  1981  for  a  comprehensive 
compilation).  Most  laboratory  data  suggest  that,  at  least  for  dry  rocks,  Q  is  independent  of  fre¬ 
quency  (Birch  and  Bancroft,  1938;  Peselnick  and  Outerbridge,  1961;  Klima  et  al.,  1964;  KnopofF, 
1964;  Pandit  and  Savage,  1973;  Toksoz  et  al.,  1979;  Nur  and  Winkler,  1980;  Johnston  and  Toksoz, 
1980;  Tittman  et  al.,  1981).  Water  saturation  generally  decreases  Q  values  of  both  P  and  S  waves, 
although  the  decrease  is  much  greater  for  S- waves  than  for  P-waves. 

In  the  laboratory  Q  increases  with  increasing  confining  pressure.  However,  the  laboratory  Q 
values  at  pressures  of  2  kilobars  or  more  in  crystalline  rocks  are  generally  less  than  1000  (Klima  et 
al.,  1964;  Bradley  and  Fort,  1966;  Mason  et  al.,  1970).  It  is  only  in  the  case  of  totally  outgassed  and 
volatile  free  rocks  that  Q  values  of  2000  or  more  have  been  obtained  (Clark  et  al.,  1980;  Tittman  et 
al.,  1974). These  values  have  been  observed  in  the  completely  dry  environment  of  the  moon  (Dainty 
et  al.,  1976).  The  Earth’s  crust  is  not  free  of  water  and  volatiles  and  the  high  Q  values  cannot  be 
attributed  to  dehydration.  Thus  the  high  Q  values  need  to  be  explained. 
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Although  Q  is  independent  of  frequency  in  dry  rocks,  it  may  be  frequency  dependent  in  saturated 
and  partially  saturated  rocks  (Gardner  et  al.,  1964;  Winkler  and  Nur,  1979;  Spencer,  1981;  Tittman 
et  al.,  1981;  Toksoz  et  al,  1979).  The  saturation  may  produce  relaxation  peaks  at  certain  frequencies 
and  increase  and  decrease  of  Q  on  two  sides  of  a  peak.  A  question  we  wish  to  investigate  is  whether 
such  relaxation  phenomena  and  fluid  motions  can  explain  the  frequency  dependence  of  crustal  Q 
values  measured  from  coda  waves. 

Interpretation  of  Strong  Motion  Data 

Scattering  of  elastic  waves  propagating  in  a  heterogeneous  medium  contributes  to  the  attenuation 
of  these  waves.  Scattering  attenuation  is  not  an  energy  dissipation  mechanism,  but  only  an  energy 
redistribution  in  space  and  time,  therefore,  it  is  a  geometric  effect.  Under  the  single  scattering 
approximation,  the  scattering  attenuation  cannot  be  separated  from  the  intrinsic  attenuation.  In 
order  to  separate  these  two  attenuation  mechanisms,  we  need  to  use  the  multiple  scattering  theory. 
There  is  no  general  solution  for  the  multiple  scattering  theory.  However,  several  special  cases 
have  been  studied  (O’Doherty  and  Anstey,  1971;  Kopnichev,  1977;  Dainty  and  Toksoz,  1977,  1981; 
Richards  and  Menke,  1983;  and  Gao  et  al.,  1983a,  b).  Wu  (1984,  1985)  formulated  the  multiple 
scattering  problem  in  the  frequency  domain  using  radiative  transfer  theory.  In  the  case  of  isotropic 
scattering  with  a  point  source  in  an  infinite  random  medium,  an  exact  solution  can  be  obtained 
(Appendix  A). 

Figure  1  shows  the  distribution  of  seismic  wave  energy  with  distance  calculated  by  the  theory. 
In  the  figure,  the  distance  is  normalized  by  the  extinction  length  Le,  which  is  the  reciprocal  of  the 
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extinction  coefficient  r)e: 


Le  =  l/»)e 

He  =  la  +  rjs 

where  r)a  is  the  energy  absorption  coefficient  due  to  anelasticity  of  the  medium  and  r),  is  the 
scattering  coefficient  which  is  defined  as  the  total  scattered  power  by  a  unit  volume  of  random 
medium  per  unit  incident  power  flux  density.  Note  that  r)a  is  related  to  the  anelastic  attenuation 
coefficient  given  in  equation  (3)  by  rja  =  2 a.  In  Figure  1  the  curve  shapes  change  depending  on  the 
seismic  albedo  Bo  of  the  medium,  which  is  defined  as: 


*)> 

1*  +  *7a 


(6) 


For  the  case  of  large  albedo  (Bo  >  0.5),  i.e.  when  the  medium  is  strongly  heterogeneous,  and 
scattering  is  significant,  the  curves  are  of  arch  shape.  The  maxima  of  the  curves  depend  on  the 
extinction  coefficient  ( rje  =  r/t  +  r/a).  Therefore  it  is  possible  to  obtain  B o  and  rje  from  the  energy 
density-distance  curves,  and  thus  separate  the  scattering  effect  from  the  intrinsic  attenuation. 

The  theory  has  been  applied  to  local  earthquakes  in  Hindu  Kush  region  (Wu,  1984;  Wu  and 
Aki,  1987)  with  the  conclusion  that  the  scattering  attenuation  in  that  region  is  not  the  dominant 
factor  (Bo  <  0.5).  In  this  study,  we  look  at  the  attenuation  data  in  the  eastern  United  States  where 
anelastic  attenuation  may  be  low. 

Figures  2a  and  2b  are  the  strong  motion  data  (pseudo  velocity)  in  Northeastern  America  for 
the  case  of  /  =  5  Hz  and  1  Hz  respectively  (with  5%  damping).  The  solid  lines  in  the  figures  are 
the  best  fits  to  the  data.  If  we  assume  that  the  received  strong  motions  are  composed  of  both  the 
direct  arrivals  and  the  scattered  waves,  then  we  can  compare  curves  given  in  Figure  2  with  the 
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Table  1:  Medium  parameters  at  f—\  and  /= 5  Hz  based  on  multiple  scattering  theory 


Parameter 

f  =1  Hz 

11 

cn 

X 

N 

Lt 

15  km 

25  km 

Bo 

0.9 

0.8 

0.06  km-1 

0.03  km-1 

M= 

16.7  km 

31.3  km 

Va 

0.0067  km-1 

0.0081  km'1 

M=  !/»?<*) 

150  km 

125  km 

Q,(=  kLt) 

30 

280 

Qa{=  kLa) 

270 

1120 

•data  to  obtain  the  seismic  albedo  Bo  and  the  intrinsic  quality  factor  Qa.  In  Figure  3,  the  PSV 
data  are  corrected  for  the  geometric  spreading  (l/r  for  body  waves)  and  then  squared  to  compare 
with  the  theoretical  predictions.  The  best  theoretical  curves  are  also  drawn  in  the  Figure.  We  can 
see  that  in  the  first  100  km  the  fit  between  theory  and  data  is  generally  good  except  for  a  few 
points  which  are  very  close  to  the  source.  For  greater  distances,  the  data  gradually  deviate  from 
the  theory  and  become  flatter.  This  may  be  due  to  the  dominance  of  Lg  waves  at  great  distances. 
The  discrepancy  of  data  and  theory  at  very  close  distances  is  probably  due  to  the  non-linear  effects. 
From  these  comparisons  of  data  with  theory,  we  obtain  the  average  seismic  albedo  Bo  =  0.9  and 
the  extinction  length  Le  =  15  km  for  1  Hz  waves  and  Bo  =  0.8,  Le  =  25  km  for  5  Hz  waves.  The 
medium  parameters  based  on  these  values  are  listed  in  Table  1. 

In  Figure  4  we  plot  the  theoretical  curves  of  the  PSV-distance  relation  for  different  seismic 


albedo  Bq  when  the  extinction  length  is  fixed  at  15  km.  A  smaller  albedo  means  a  smaller  intrinsic 


Q  and  therefore  has  a  steep  decrease  of  amplitude  with  distance.  Figure  5  shows  different  curves  of 
different  albedos  when  the  intrinsic  Q  is  fixed  at  Qa  =  1350.  We  can  see  that  the  strong  scattering 
will  make  the  apparent  attenuation  much  bigger  than  the  intrinsic  attenuation  when  the  distance  is 
larger  than  the  absorption  extinction  length,  La.  However,  the  amplitude  change  is  not  exponential 
for  small  distances. 

Results  given  in  Table  1  give  a  good  fit  to  the  data  with  a  consistent  set  of  parameters  at  /  = 
1  Hz  and  /  =  5  Hz.  They  suggest  a  frequency  dependent  anelastic  Qa  =  270/°  88.  As  discussed  in 
the  first  section  while  reviewing  the  laboratory  data,  such  variation  of  Q  with  frequency  cannot  be 
explained  without  a  relaxation  mechanism.  In  the  crust,  the  fluids  may  provide  such  a  mechanism. 

Attenuation  of  Rg 

The  data  set  used  consists  of  recordings  of  blasts  for  a  refraction  experiment  in  southeastern  Maine 
conducted  by  the  USGS  (Doll  et  al.,  1986).  The  largest  phase  from  these  near  surface  sources, 
and  hence  the  phase  with  the  peak  amplitude,  is  Rg,  the  fundamental  mode  Rayleigh  wave.  Some 
examples  of  the  seismograms  are  shown  in  Figure  6.  The  attenuation  of  Rg  can  provide  useful 
information  first,  because  it  shows  the  affect  of  the  attenuation  mechanisms  on  an  isolated,  relatively 
short  pulse,  and  second,  it  gives  information  about  the  depth  dependence  of  Q  from  the  frequency 
dependence.  The  frequencies  observed  on  the  records  for  Rg  are  1-10  Hz.  To  calculate  Q,  equation 
(1)  is  rewritten  for  Rg  in  terms  of  the  attenuation  coefficient  7(0/)  as 

ln(v/rl/(r,w))  =  lnS(w)  -  r  •  7(w)  (7) 
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A  windowed  portion  of  seismograms  at  different  ranges  r  from  a  single  shot  is  Fourier  transformed 
to  find  U(r,u ).  Since  the  spectrum  window  contains  the  entire  Rg  wave  train,  there  is  no  correction 
for  dispersion.  The  source  term  S(u>)  is  common  to  all  seismograms.  Accordingly,  the  slope  of  a 
least  square  fit  straight  line  to  the  left  side  of  (7)  will  determine  7.  The  total  Qt  in  the  sense  of 
equation  (2)  may  be  calculated  from  Qt  =  w/(2^v).  Note  that  we  distinguish  the  experimentally 
determined  attenuation  coefficient  7' from  the  anelastic  attenuation  coefficient  a  in  equation  (3), 
since  the  two  are  not  in  general  equal  by  equation  (2). 

In  Figure  7  we  show  average  Q  for  the  area  together  with  models  and  fits  to  the  data  under  three 
different  assumptions  concerning  the  frequency  dependence  of  Q.  Discussing  the  measured  values 
first,  all  Q  values  are  less  than  50,  and  this  indicates  strong  attenuation  in  the  uppermost  crust. 
The  Rg  Q  values  are  substantially  lower  than  those  derived  by  Pulli  (1984)  from  coda  discussed 
below  and  have  a  different  frequency  dependence:  Pulli  obtained  a  Q  of  460  at  1  Hz,  increasing 
with  frequency.  It  should  be  noted  that  at  all  frequencies  the  coda  results  are  probably  averages  for 
the  whole  crust  while  Rg  Q  is  sensitive  to  near  surface  Q  and  has  a  depth  dependence  that  varies 
with  frequency.  While  it  is  not  possible  to  separate  Qt  into  Qa  and  Q the  Qt  values  are  similar  to 
Q,  in  the  previous  section  and  suggest  that  scattering  is  important,  although  Qa  is  probably  also 
low,  as  discussed  below. 

Rg  Q  at  any  frequency  is  a  weighted  average  over  depth  of  both  compressional  and  shear  Q,  but 
mainly  shear  Q.  This  weighting  depends  on  frequency,  since  the  lower  the  frequency  the  greater 
the  penetration  of  the  Rg  motion.  It  then  becomes  difficult  to  separate  frequency  variation  of  Rg  Q 
due  to  depth  dependence  from  frequency  dependence  at  a  given  depth.  As  a  first  step  in  attacking 
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this  difficulty,  Figure  7  gives  results  for  three  different  cases  in  which  the  frequency  dependence 
of  shear  Q  is  fixed  at  all  depths  as  Qk  =  Qokfn >  where  k  is  the  layer  index  of  a  layered  model 
and  n  =  0,0.5, 1  respectively,  spanning  the  observed  range  (Toksoz  et  al.,  1988).  In  Figure  7  Qok 
( i.e. ,  Q  at  1  Hz)  is  plotted  as  a  function  of  depth  with  the  fit  to  the  observations  shown  above. 
The  three  choices  of  n  can  be  thought  of  as  representing  attenuation  due  to  microcracks  and  grain 
boundary  sliding  (n  =  0;  Toksoz  and  Johnston,  1981),  attenuation  due  to  fluid  flow  in  major  open 
fractures  (n  =  0.5;  Toksoz  et  al.,  1987),  and  attenuation  due  to  scattering  (n  =  1;  Dainty,  1981). 
Other  mechanisms  of  attenuation  are  also  possible,  but  the  three  presented  here  must  certainly  be 
examined  in  the  low  pressure,  highly  heterogeneous  environment  of  the  upper  crust. 

For  all  of  the  models  Q  may  only  be  determined  down  to  a  depth  of  about  3  km  due  to  the 
limited  penetration  of  Rg.  From  Figure  7,  the  constant  Q  model  (n  =  0)  has  lower  Q  at  depth 
than  near  the  surface.  This  contradicts  laboratory  data  and  investigations  in  other  regions  and  we 
reject  this  model.  For  the  case  n  =  1  data  at  high  frequencies  cannot  be  fit  because  of  the  strong 
frequency  dependence  of  Q,  although  a  reasonable  fit  is  found  below  2  Hz.  A  better  fit  is  obtained 
for  n  =  0.5,  except  below  1.5  Hz. 

To  extend  the  attenuation  profile  deeper  into  the  crust  we  use  Rayleigh  wave  data  from  central 
North  America  in  the  period  range  4  to  20  sec  (Herrmann  and  Mitchell,  1975;  Mitchell,  1980).  A 
Q  structure  that  fits  both  data  sets  is  shown  in  Figure  8.  In  the  top  2  km  of  this  model  the  Q 
structure  is  determined  by  the  high  frequency  (/  >  0,85  Hz)  data  presented  above;  the  frequency 
dependence  is  Q  =  Qo/'j5-  The  deeper  Q  structure  is  controlled  by  the  longer  period  data  of 
Mitchell  (1980).  To  estimate  the  frequency  dependence,  we  note  that  the  attenuation  (Q_1)  due  to 
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scattering  or  fluid  flow  is  peaked  at  a  critical  frequency.  Below  the  critical  frequ<  :y  attenuation 
decreases  rapidly  with  decreasing  frequency:  Q-1  is  proportional  to  /  for  fluid  flow  (Toksoz  et  al., 
1987)  and  to  /2  for  scattering  (Frankel  and  Clayton,  1986).  The  effect  is  make  frequency  dependent 
Q  negligably  small  at  low  frequency  and  intrinsic  anelasticity  dominates.  We  assume  Q  ~  Q o  for 
/  <  0.85  Hz  in  the  inversion  of  the  combined  data  sets. 

The  final  attenuation  model  has  low  Q  at  shallow  depth  and  high  Q  (l  1000)  below  about  5 
km.  Constraining  the  details  of  the  Q  profile  in  the  transition  zone  requires  Rg  attenuation  data 
for  0.2  <  /  <  1.0  Hz. 

Coda  Q 

Since  the  introduction  of  the  method  by  Aki  and  Chouet  (1975),  measurement  of  Q  from  the  decay 
of  coda  has  become  a  standard  technique.  There  has,  however,  been  some  question  regarding  what 
has  been  measured.  To  help  answer  this  question  we  have  compared  our  results  for  strong  motion 
Q  to  the  work  of  Pulli  (1984)  on  coda  Q  in  New  England.  Pulli  noted  that  there  was  an  apparent 
change  in  the  decay  rate  of  the  coda  at  about  100  s  after  earthquake  origin  time,  leading  to  different 
estimates  of  coda  Q,  Qc.  Before  100  s 

Qc  =  140/095  (8) 

where  /  is  frequency,  while  at  times  longer  than  100  s 

Qc  =  660 f04  (9) 
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The  average,  ignoring  the  decay  rate  change,  is 


Qc  =  460/° 4  (10) 

Note  that  for  all  cases  Qc  is  similar  to  Qa  previously  obtained  from  the  strong  motion  data  and 
greater  than  Q, .  The  Qc  found  for  shorter  times  unquestionably  refers  to  the  crust  and  is  closest 
to  the  Qa  results  for  the  strong  motion  data,  but  may  still  have  some  contribution  from  scattering 
(Dainty  et  al.,  1987).  The  longer  time  results  are  consistent  with  the  investigation  of  Singh  and 
Herrman  (1983)  covering  the  continental  U.  S.,  although  it  is  not  clear  whether  they  refer  to  the 
crust  only  or  to  the  whole  lithosphere.  It  is  less  likely  that  scattering  is  included  in  Qc  at  long 
times.  Of  particular  interest  are  the  frequency  exponents  n  =  0.95  and  n  =  0.4.  Values  close  to  0.5 
have  often  been  obtained  in  other  studies  if  Q  at  1  Hz  (Qo)  is*  not  too  high.  On  the  other  hand, 
if  Qo  is  high,  then  the  frequency  exponent  is  often  low,  for  example,  Singh  and  Herrman  (1983) 
quote 

Qc  =  1000/° 2  (11) 

for  the  central  U.  S.  In  tectonic  regions  and  at  short  coda  times  Qo  is  often  smaller  (~  100) 
and  n  may  range  up  to  1  (Singh  and  Herrman,  1983;  Herraiz  and  Espinosa,  1987;  Dainty  et  al., 
1987).  In  this  paper  we  will  make  the  preliminary  interpretation  that  coda  Q  in  stable  regions  is 
representative  of  anelastic  Q  because  of  the  similarity  with  the  strong  motion  Qa  noted,  and  that 
because  coda  Q  can  be  measured  easily  and  reliably  it  gives  the  best  estimate  of  the  frequency 
dependence  of  the  aneKstic  Q. 

One  reason  for  assuming  coda  Q  is  a  measure  of  average  crustal  Q,  at  least  in  some  cases,  is 
the  similarity  of  coda  Q  and  Lg  Q.  Recently  Campillo  et  al.  (1985)  have  demonstrated  this  in  a 
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detailed  study  in  central  France  of  both  Lg  Q  and  Q  in  the  early  coda,  obtaining  for  both  cases 

Q  =  300/°  5  (12) 

This  is  not  very  different  from  the  average  result  of  Pulli  (1984). 

Model  of  Crustal  Scattering  and  Attenuation 

Our  studies  of  scattering  and  anelastic  attenuation  with  distance  at  frequencies  of  1-10  Hz  indicate 
that  the  relative  effect  of  scattering  is  strong.  The  result  of  this  is  a  complex  interplay  between  the 
type  of  measurement  made  and  the  appropriate  attenuation  relation.  Measurements  of  quantities 
such  as  pulse  amplitude  will  feel  the  full  force  of  scattering  attenuation  and  will  decay  rapidly  with 
Q  dominated  by  scattering.  Measurements  3uch  as  coda  will  decay  less  rapidly  because  some  of 
the  scattered  energy  is  rescattered  back  into  the  measurement  window.  A  crucial  parameter  in 
the  description  of  these  effects  is  the  albedo,  which  is  the  ratio  of  attenuation  due  to  scattering 
to  total  attenuation.  In  New  England  we  found  an  albedo  of  0.8-0.9,  indicating  the  importance 
of  scattering.  This  strong  scattering  means  that  coda  Q  is  probably  most  closely  identified  with 
anelastic  crustal  Q,  which  appears  to  be  frequency  dependent.  The  results  from  Rg  indicate  that 
there  is  a  low  Q  surface  layer  about  2  km  thick  with  Q  <  100. 

Laboratory  data  and  observations  made  in  boreholes  help  to  achieve  an  understanding  of  the 
physical  mechanisms  contributing  to  the  attenuation  of  seismic  waves  at  regional  distances.  Labo¬ 
ratory  measurements  of  attenuation  in  “wet"  rocks  give  relatively  low  intrinsic  Q  (Q  <  200)  even  at 
pressures  of  2  kb,  and  do  not  indicate  that  Q  values  would  increase  significantly  with  pressure.  A 
primary  mechanism  for  this  attenuation  is  friction  along  grain  boundaries  and  microcracks  (Walsh, 
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1966;  Toksoz  and  Johnston,  1981).  Recent  studies  indicate  that  in  “wet”  rocks  microcracks  heal  at 
pressures  of  about  2  kb  and  temperatures  of  400°C  (Smith  and  Evans,  1984;  Rovetta  et  al.,  1987). 
Thus,  microcracks  would  not  contribute  to  attenuation  in  the  crust  below  a  few  kilometers  and 
intrinsic  Q  below  this  depth  is  likely  to  be  high  ( Q  >  1000). 

Data  from  wells  drilled  in  crystalline  rocks  shows  the  presence  of  both  major  and  minor  fractures. 
The  fluid  flow  in  such  fractures  can  contribute  significantly  to  attenuation  and  give  a  frequency 
dependent  Q  (Toksoz  et  al.,  1987).  Above  a  critical  frequency,  which  could  be  as  low  as  1  or  2  Hz,  Q 
is  proportional  to  the  square  root  of  frequency  (/°  5).  Below  the  critical  frequency  Q  is  proportional 
to  inverse  frequency  (/-1).  There  is  an  even  stronger  frequency  dependence  of  attenuation  due  to 
scattering  with  a  peak  between  /  =  aj A  and  /  =  2o/A,  where  a  is  the  correlation  length  of  random 
scatterers  and  A  is  the  wavelength  of  seismic  waves.  Above  and  below  this  peak  Q  is  proportional  to 
frequency  (/*)  and  inverse  frequency  squared  (/~2)  respectively  (Wu,  1982;  Dainty,  1984;  Frankel 
and  Clayton,  1986). 

With  all  these  mechanisms  included,  a  conceptual  attenuation  and  scattering  model  for  the 
crust  can  be  summarized.  In  the  shallow  crust  friction  mechanisms  (constant  Q),  fluid  flow  and 
scattering  all  contribute  to  attenuation.  With  increasing  depth  the  effects  of  friction  and  fluid  flow 
decrease  rapidly.  Recent  results  of  Flatte  and  Wu  (1988)  indicate  that  the  strong  scattering  regime 
may  extend  to  10-15  km  depth. 

In  terms  of  numbers  applicable  to  average  continental  crust  in  non-tectonic  areas  without  a 
thick  section  of  recent  sediments,  the  following  approximate  model  can  be  used: 

1.  At  shallow  depth  (h  <  5  km)  low  intrinsic  (anelastic)  Q  with  frequency  dependence  /Oo.  A 
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typical  value  would  be  Q  =  100/° 5 .  Scattering  is  also  important  in  this  region. 

2.  In  the  lower  crust  high  anelastic  Q  (Qo  >  1000).  The  influence  of  scattering  may  be  moderate 
and  remains  to  be  determined  in  this  region. 

These  values  are  consistent  with  models  proposed  by  Campillo  et  al.  (1985)  for  frequency  dependent 
attenuation  in  the  crust  in  central  France  based  on  Lg  waves,  and  with  the  results  shown  in  Figure  8. 
This  model  is  also  consistent  with  the  model  of  lithospheric  heterogeneities  derived  by  Flatte  and  Wu 
(1988)  using  travel  time  and  amplitude  fluctuations  at  the  Norwegian  Seismic  Array  (NORSAR). 
In  their  model  the  upper  crust  to  a  depth  of  15  km  is  rich  in  small  (less  than  1  km  scale  size) 
heterogeneities,  while  below  that  both  small  and  intermediate  scale  heterogeneities  are  found  to 
about  200  km  depth.  The  strongly  heterogeneous  nature  of  the  upper  crust  has  also  been  confirmed 
by  reflection  and  -efraction  profiling  (Mooney  and  Brocher,  1987). 

Acknowledgements 

This  work  was  supported  by  the  United  States  Geological  Survey  under  contract  number  14-08- 
0001-G1092  and  by  the  Advanced  Research  Project  Agency  of  the  Department  of  Defense  and 
monitored  by  the  Air  Force  Geophysical  Laboratory  under  contract  number  F19628-86-K-0004. 
The  views  expressed  in  this  report,  however,  are  solely  those  of  the  authors  and  do  not  necessarily 
represent  the  views  of  the  United  States  Geological  Survey,  the  Advanced  Research  Projects  Agency, 
the  Air  Force  Geophysical  Laboratory,  or  the  United  States  Government. 


62 


REFERENCES 


Aki,  K.  (1980),  Attenuation  of  shear  waves  in  the  lithosphere  for  frequencies  from  0.05  to  25  Hz, 
Phys.  Earth  Planet.  Inter.  21,  50-60. 

Aki,  K.,  and  Chouet,  B.  (1975),  Origin  of  coda  waves:  source,  attenuation  and  scattering  effects,  J. 
Geophys.  Res.  80,  3322-3342. 

Aki,  K.,  and  Richards,  P.G.,  Quantitative  seismology.  Theory  and  methods  (W.H.  Freeman  and 
Co.,  San  Francisco  1980). 

Birch,  F.,  and  Bancroft,  D.  (1938),  Elasticity  and  internal  friction  in  a  long  column  of  granite,  Bull. 
Seism.  Soc.  Am.  28,  243-254. 

Bradley,  J.J.,  and  Fort,  A.N.,  Jr.,  Internal  friction  in  rocks,  In  Handbook  of  Physical  constants  (ed. 
Clark,  S.P.,  Jr.)  (1966)  pp.  175-193. 

Campillo,  M.,  Plantet,  J.-L.,  and  Bouchon,  M.  (1985),  Frequency- dependent  attenuation  in  the  crust 
beneath  central  France  from  Lg  waves:  data  analysis  and  numerical  modeling,  Bull.  Seism.  Soc. 
Am.  75,  1395-1411. 

Clark,  V.A.,  Spencer,  T.W.,  Tittmaun,  B.R.,  Ahlberg,  L.A.,  and  Coombe,  L.T.  (1980),  Effect  of 
volatiles  on  attenuation  Q  _1  and  velocity  in  sedimentary  rocks,  J.  Geophys.  Res.  85,  5190-5198. 
Dainty,  A.M.  (1981),  A  scattering  model  to  explain  seismic  Q  observations  in  the  lithosphere  between 
1  and  80  Hz,  Geophys.  Res.  Lett.  8,  1126-1128. 

Dainty,  A  M.  (1984),  High-frequency  acoustic  backscattering  and  seismic  attenuation,  J  Geophys. 
Res.  89,  3172-3176. 

Dainty,  AM.,  and  Toksoz.  M  N  (1977).  Elastic  wave  propagation  in  a  highly  scattering  medium.  .4 


63 


diffusion  approach,  J.  Geophys.  I>S,  375-388. 

Dainty,  AM.,  and  Toksoz,  M.N.  (1981),  Seismic  codas  on  the  Earth  and  the  Moon:  a  comparison, 
Phys.  Earth  and  Plan.  Int.  26 ,  250-260. 

Dainty,  A.M.,  Toksoz,  M.N.,  and  Stein,  S.  (1976),  Seismic  investigation  of  the  lunar  interior,  Lunar 
Science  Conf.,  7th  Proc.,  Geochim.  et  Cosmochim.  Acta,  suppl.  7  5,  3057-3075. 

Dainty,  A.M.,  Duckworth,  R.M.,  and  Tie,  A.  (1987),  Attenuation  and  backscattering  from  local  coda, 
Bull.  Seism.  Soc.  Am.  77,  1728-1747. 

Doll,  W.E.,  Luetgert,  J.H.,  and  Murphy,  J.M.  (1986),  Seismic  refraction  and  wide-angle  reflection 
in  the  northern  Appalachians;  the  Chain  Lakes  Massif,  northwest  Maine  (abst.),  EOS  67,  312. 

Flatte,  S.M.,  and  Wu,  R.S.  (1988),  Small-scale  structure  in  the  lithosphere  and  asthenosphere  de¬ 
duced  from  arrival  time  and  amplitude  fluctuations  at  NORSAR ,  J.  Geophys.  Res.,  in  press. 

Frankel,  A.,  and  Clayton,  R.W.  (1986),  Finite  difference  simulations  of  seismic  scattering:  im¬ 
plications  for  the  scattering  of  short-period  seismic  waves  in  the  crust  and  models  of  crustal 
heterogeneity ,  J.  Geophys.  Res.  91,  6465-6489. 

Gao,  L.S.,  Lee,  L.C..  Biswas,  N.H.,  and  Aki,  K.  (1983a),  Comparison  of  the  effects  between  single 
and  multiple  scattering  on  coda  waves  for  local  earthquakes,  Bull  Seism.  Soc.  Am.  73,  373-389. 

Gao,  L.S.,  Lee,  L.C..  Biswas,  N.H.,  and  Aki,  K.  (1983b),  Effects  of  multiple  scattering  on  coda  waves 
in  three-dimensional  medium,  PAGEOPH  121,  3  -15. 

Gardner,  G.H.F.,  Wyllie,  M.R.J.,  and  Droschack,  D.M.  (1961),  Ey  cts  of  pressure  and  fluid  satu¬ 
ration  on  the  attenuation  of  elastic  waves  in  sands,  J.  Petr.  Tech.,  189-198. 

Gupta,  I.N.,  Burnetti,  A.,  McElfresh,  T.W.,  Von  Seggern,  D.H.,  and  Wagner,  R.A.  (1983),  Lateral 


64 


variations  in  attenuation  of  ground  motion  in  the  eastern  United  States  based  on  propagation 
of  Lg,  NUREG/CR-3555,  U.S.  Nuclear  Regulatory  Commission,  Washington,  D.C. 

Herraiz,  M.,  and  Espinosa,  A.  F.  (1987),  Coda  waves:  a  review,  PAGEOPH  125,  499-577. 

Herrmann,  R.B.  (1980),  Q  estimates  using  the  coda  of  local  earthquakes,  Bull.  Seism.  Soc.  Am.  70, 
447-468. 

Herrmann,  R.B.,  and  Mitchell,  B.J.  (1975),  Statistical  analysis  and  interpretation  of  surface-wave 
anelastic  attenuation  data  for  the  stable  interior  of  North  America,  Bull.  Seism.  Soc.  Am.  65, 
1115-1128. 

Johnston,  D.H.,  and  Toksoz,  M.N.  (1980),  Ultrasonic  P  and  S  wave  attenuation  in  dry  and  saturated 
rocks  under  pressure,  J.  Geophys.  Res.  85,  925-936. 

Klima,  K.,  Vanek,  J.,  and  Pros,  Z.  (1964),  The  attenuation  of  longitudinal  waves  in  diabase  and 
greywacke  under  pressure  up  to  4  kilobars,  Studia  Geoph.  et  Geod.  8,  247-254. 

Knopoff,  L.  (1964),  Q,  Rev.  Geophys.  2,  625-660. 

Kopnichev,  Y.F.  (1977),  The  role  of  multiple  scattering  in  the  formation  of  a  seismogram  tail,  Izv. 
Acad.  Sci.,  USSR,  Phys.  Solid  Earth  IS,  394-398. 

Mason,  W.P.,  Beshers,  D.N.,  and  Kuo,  J.T.  (1970),  Internal  friction  in  Westerly  granite:  relation 
to  dislocation  theory,  J.  Appl.  Phys.  41,  5206-5209. 

Mitchell,  B.J.  (1980),  Frequency  dependence  of  shear  wave  internal  friction  in  the  continental  crust 
of  eastern  North  America,  J.  Geophys.  Res.  85,  5212-5218. 

Mooney,  W.D.,  and  Brocher,  T.M.  (1987),  Coincident  seismic  reflection/refraction  studies  of  the 
continental  lithosphere:  a  global  review,  Rev.  Geophys.  25,  723-742. 


65 


Nur,  A.,  and  Winkler,  K.  (1980),  The  role  of  friction  and  fluid  flow  in  wave  attenuation  in  rocks 
(abst.),  Geophysics  ^5,  591-592. 

O’Doherty,  R.F.,  and  Anstey,  N.A.  (1971),  Reflections  on  amplitudes,  Geophys.  Prospect.  19,  430- 
458. 

Pandit,  B.I.,  and  Savage,  J.C.  (1973),  An  experimental  test  of  Lomnitz’s  theory  of  internal  friction 
in  rocks,  J.  Geophys.  Res.  78,  6097-6099. 

Peselnick,  L.,  and  Outerbridge,  W.F.  (1961),  Internal  friction  in  shear  and  shear  modulus  of  Solen- 
hofen  limestone  over  a  frequency  range  of  107  cycles  per  second,  J.  Geophys.  Res.  66,  581-588. 

Pulli,  J.  (1984),  Attenuation  of  coda  waves  in  New  England,  Bull.  Seism.  Soc.  Am.  7f,  1149-1166. 

Rautian,  T.G.,  and  Khalturin,  V  I.  (1978),  The  use  of  the  coda  for  the  determination  of  the  earth¬ 
quake  source  spectrum,  Bull.  Seism.  Soc.  Am.  68,  1809-1827. 

Richards,  P.G.,  and  Menke,  W.  (1983),  The  apparent  attenuation  of  a  scattering  medium,  Bull. 
Seism.  Soc.  Am.  73,  1005-1021. 

Roecker,  S.W.,  Tucker,  B.,  King,  J.,  and  Hatzfield,  D.  (1982),  Estimates  of  Q  in  central  Asia  as  a 
function  of  frequency  and  depth  using  the  coda  of  locally  recorded  earthquakes,  Bull.  Seism.  Soc. 
Am.  72,  129-149. 

Rovetta,  M.R.,  Blacic,  J.D.,  and  Delaney,  J.R.  (1987),  Microfracture  and  crack  healing  in  experi¬ 
mentally  deformed  peridotite,  J.  Geophys.  Res.  92,  12902-12910. 

Singh,  S.  (1985),  Lg  and  coda  wave  studies  of  Eastern  Canada,  Ph  D.  thesis,  Saint  Louis  University, 
Missouri. 

Singh,  S.,  and  Hermann,  R.B.  (1983),  Regionalization  of  crustal  coda  Q  in  the  continental  United 


66 


I 


States,  J.  Geophys.  Res.  88,  527-538. 

Smith,  D.L.,  and  Evans,  B.  (1984),  Dtffusional  crack  healing  in  quartz,  J.  Geophys.  Res  89,  4125 — 

4135. 

Spencer,  J.W.,  Jr.  (1981),  Stress  relaxations  at  low  frequencies  in  fluid  saturated  rocks:  attenuation  | 

and  modulus  dispersion,  J.  Geophys.  Res.  86,  1803-1812. 

Tittmann,  B.R.,  Iiousley,  R.M.,  Alers,  G.A.,  and  Cirlin,  E.H.  (1974),  Internal  friction  in  rocks 
and  its  relationship  to  volatiles  on  the  moon,  Lunar  Science  Conf.,  5th  Proc.,  Geochim.  et 
Cosmochim.  Acta,  suppl.  5  8,  2913-2918. 

Tittmann,  B.R.,  N'a  'ler,  H.,  Clark,  V.A.,  Ahlberg,  L.A.,  and  Spencer,  T.W.  (1981),  Frequency 
dependence  of  seismic  dissipation  in  saturated  rocks,  Geophys.  Res.  Lett.  8,  36-38. 

Toksoz,  M.N.,  and  Johnston,  D.H.,  (eds),  (1981),  Seismic  wave  attenuation ,  Geophysics  reprint 
series  No.  2,  Society  of  Exploration  Geophysicist. 

Toksoz,  M.N.,  Johnston,  D.H.,  and  Timur,  A.  (1979),  Attenuation  of  seismic  waves  in  dry  and 
saturated  rocks.  I.  Laboratory  measurements,  Geophysics  44,  681-690. 

Toksoz,  M.N.,  \\u,  R.S.,  and  Schmitt,  D.P.  (1987),  Physical  mechanisms  contributing  to  seismic 
attenuation  in  the  crust,  Proc.  NATO  ASI  ‘Strong  Ground  Motion  Seismology’,  Ankara,  Turkey, 

225-247. 

Toksoz,  M.N.,  Reiter,  E.,  and  Nlandal,  B.  (1988),  Seismic  attenuation  in  the  crust,  Proc.  10th  Ann. 
DARPA/AFGL  Res.  Symp.,  127-134. 

Walsh,  J.B.  (1966),  Seismic  wave  attenuation  in  rock  due  to  fracture,  J.  Geophys.  Res.  17,  2591- 
2599. 


67 


Winkler,  K.,  and  Nur,  A.  (1979),  Pore  fluids  and  seismic  attenuation  in  rocks,  Geophys.  Res.  Lett. 
6,  1-4. 

Wu,  R.S.  (1982),  Attenuation  of  short  period  seismic  waves  due  to  scattering,  Geophys.  Res.  Lett. 
9,  9-12. 

Wu,  R.S.  (1984),  Multiple  scattering  and  energy  transfer  of  seismic  waves  and  the  application  of  the 
theory  to  Hindu  Kush  region,  In  Seismic  wave  scattering  and  the  small  scale  inhomogeneities  in 
the  lithosphere  (Chapter  4,  Ph.D.  thesis,  Mass.  Inst.  Tech.,  Cambridge,  MA). 

Wu,  R.S.  (1985),  Multiple  scattering  and  energy  transfer  of  seismic  waves.  Separation  of  scattering 
effect  from  intrinsic  attenuation.  I.  Theoretical  modelling,  Geophys.  J.  R.  Astr.  Soc.  82,  57-80. 

Wu,  R.S.,  and  Aki,  K.  (1987),  Multiple  scattering  and  energy  transfer  of  seismic  waves — separation 
of  scattering  effect  from  intrinsic  attenuation.  II.  Application  of  the  theory  to  the  Hindu  Kush 
region.  This  issue. 


68 


APPENDIX  A 


The  elastic  wave  energy  received  by  an  isotropic  point  receiver  in  a  random  heterogeneous  medium 
can  be  represented  by  the  average  energy  density  E(rf  at  that  point.  The  energy  density  in  radiative 
transfer  theory  is  defined  as  (see  Wu,  1985): 

E{r)  =  ~  jjir^d  n  (A-  1) 

where  C  is  the  wave  velocity,  and  /(r,n)  is  the  specific  intensity  or  directional  intensity.  It  gives 
the  power  flowing  within  a  unit  solid  angle  in  the  direction  Cl  (Cl  is  the  unit  vector)  received  by 
a  unit  area  perpendicular  to  f2,  in  a  unit  frequency  band.  The  specific  intensity  is  defined  for  a 
frequency  w,  which  is  omitted  in  the  notation.  Since  the  P  wave  energy  is  much  smaller  than  the 
S  wave  energy  for  earthquakes,  we  consider  here  /(r,  f2)  as  only  the  S  wave  energy  by  neglecting 
the  mode  converted  energy  from  P  waves.  We  assume  here  also  that  the  wave  energy  described  by 
/(r,  f2)  is  depolarized,  i.e.  the  energy  is  equally  partitioned  between  the  two  orthogonal  components 
of  S  waves.  This  agrees  generally  with  the  observations. 

From  the  radiative  transfer  equations  we  can  obtain  the  equation  for  the  transfer  energy  density: 

E(r)  =  Eine ~r,'t  -  Jv  [r7,E(nJ  *  e(n,  H)]  G0(r  -  rJ)dVl  (A  -  2) 

where  E{n  is  the  incident  field  and, 

e(H.n)=^W(rLtn)  (A -3) 

is  the  source  energy  density  function,  where  W  is  the  source  power  spectral  density,  and 

e-n-R 

AnR“  4 7r  r  -  r,  2 


(A  -  4) 


G  (n  -  L) 
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and 

/(*,B0)=|[l-^tanh-1(I)]‘-r^j  j  (A -10) 

The  first  term  in  Equation  (A-7)  is  the  diffuse  term  Ej  and  the  second  term  is  the  coherent  term 

Ec. 

Note  that  the  diffuse  multiplier  do  is  always  less  tha..  1.  When  distance  r  is  large,  especially 
For  large  Bo,  the  diffuse  term  becomes  dominant,  and  E(r)  will  be  approximately  an  exponential 
decay  with  an  apparent  attenuation  coefficient  d^Ve,  which  is  less  than  the  extinction  coefficient  17,. 
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The  degree  of  reduction  depends  on  the  albedo  Bq.  Figure  1  shows  the  energy  density  distribution 
with  distances  for  different  albedo  values. 
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Figure  Captions 


Figure  1.  Normalized  energy  distribution  curves  corrected  for  spherical  spreading,  4ir r2 E(r)  as  a 
function  of  normalized  distance  De  =  r/Le  where  Le  is  the  extinction  length  defined  by  Equation 
6  in  the  text. 

Figures  2a, b.  Ground  velocity  (PSV)  at  5  Hz  (3a)  and  1  Hz  (2a)  as  a  function  of  distance  for  events 
in  northeastern  United  States  and  Eastern  Canada.  Values  normalized  to  a  common  magnitude. 
Data  are  from  compilation  of  Risk  Engineering,  Inc.,  under  EPRI  sponsorship.  The  solid  line 
in  each  case  is  a  best  fit  to  data. 

A  11-1-82,  New  Brunswick,  A {  -  b  =  5.5,  ECTN  data 

O  19-1-82,  New  Hampshire,  Afj  =  4.8,  strong  motion  data  and  ECTN 

□  31-3-82,  New  Brunswick,  Mj  =  4.8.  strong  motion  data  and  ECTN 

O  6-5-82,  New  Brunswick,  Afj,  =  4.0,  strong  motion  data 

▲  16-6-82,  New  Brunswick,  =  4.6,  strong  motion  data  and  ECTN 

0  7-10-83,  Adirondacks,  New  York,  Mf,  —  5.6,  ECTN 

■  11-10-83,  Ottawa,  Canada,  Aft,  =  4.1,  ECTN 

Figure  3a, b.  Match  between  the  multiple  scattering  model  and  the  observed  ground  motion  data  as 
a  function  of  radial  (epicentral)  distance  R,  at  frequencies  5  Hz  (3a)  and  1  Hz  (3b).  PSV  curves 
ere  the  best  fit  curves  of  Figures  2a, b.  ( PSV  ■  R  10)  and  (PSV  ■  R/10)2  are  calculated  from 

PSV  curves.  Model  and  data  curves  are  fit  in  the  distance  range  of  R  -  10  to  100  km  where 

model  approximations  are  valid.  Results  of  the  fitting  are  given  in  Table  1. 
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Figure  4.  Sensitivity  of  theoretical  curves  (Power  versus  radial  distance)  at  /  =  5  Hz,  to  different 
model  parameters. 

Figure  5.  Sensitivity  of  theoretical  curves  to  albedo  ( Ba )  values  at  /  =  5  Hz  as  a  function  of  distance. 
Fixed  parameters  are  Le  =  15  km,  Qa  —  1350. 

Figure  6.  Typical  seismograms  used  in  the  Rg  study.  Source-receiver  distances  are  indicated.  Rg  is 
the  large  low-frequency  phase. 

Figure  7.  Observed  Rg  Q  versus  frequency  and  fit  (top)  and  shear  Q  structures  used  to  fit  the  data; 
Qok>  i-e->  Q  at  1  Hz  is  shown  (bottom).  Left:  frequency  independent  Q  —  Qo*  for  each  layer. 
Center:  /°  5  dependence,  Q  =  Qokf 0  5  for  each  layer.  Right:  / 1  dependence,  Q  =  Qokf 1  f°r  eac^ 
layer. 

Figure  8.  Q  model  for  the  crust  derived  from  short  period  Rg  data  of  this  study  and  longer  period 
data  of  Mitchell  (1980).  Top:  fit  to  data.  Bottom:  Q  model.  Top  2.5  km  of  the  model  is  the 
same  as  Figure  7  (center)  with  frequency  dependence  of  Q  =  Qofn,  where  n  =  0.5  for  /  >  0.85 
Hz,  and  n  =  0  for  /  <  0.85  Hz.  Below  2.5  km  Q  is  independent  of  frequency. 
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